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);

FileIO/AsciiRasterInterface.cpp

+/**
+ * @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

FileIO/AsciiRasterInterface.h

+/**
+ * @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_ */

FileIO/CMakeLists.txt

 # Source files
 # GET_SOURCE_FILES(SOURCES_FILEIO)
 SET( SOURCES
+	AsciiRasterInterface.h
+	AsciiRasterInterface.cpp
 	GMSInterface.h
 	GMSInterface.cpp
 	GMSHInterface.h

FileIO/GMSHInterface.cpp

@@ -103,130 +103,76 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 {
 	std::string line;
 	std::ifstream in(fname.c_str(), std::ios::in);
-	getline(in, line); // Node keyword
+	if (!in.is_open())
+	{
+		WARN ("GMSHInterface::readGMSHMesh() - Could not open file %s.", fname.c_str());
+		return nullptr;
+	}
+
+	getline(in, line); // $MeshFormat keyword
+	if (line.find("$MeshFormat") == std::string::npos)
+	{
+		in.close();
+		WARN ("No GMSH file format recognized.");
+		return nullptr;
+	}
+
+	getline(in, line); // version-number file-type data-size
+	if (line.substr(0,3).compare("2.2") != 0) {
+		WARN("Wrong gmsh file format version.");
+		return nullptr;
+	}
+
+	if (line.substr(4,1).compare("0") != 0) {
+		WARN("Currently reading gmsh binary file type is not supported.");
+		return nullptr;
+	}
+	getline(in, line); //$EndMeshFormat
+
 	std::vector<MeshLib::Node*> nodes;
 	std::vector<MeshLib::Element*> elements;
-
-	if (line.find("$MeshFormat") != std::string::npos)
+	std::map<unsigned, unsigned> id_map;
+	while (line.find("$EndElements") == std::string::npos)
 	{
-		getline(in, line); // version-number file-type data-size
-		if (line.substr(0,3).compare("2.2") != 0) {
-			WARN("Wrong gmsh file format version.");
-			return nullptr;
-		}
-		if (line.substr(4,1).compare("0") != 0) {
-			WARN("Currently reading gmsh binary file type is not supported.");
-			return nullptr;
-		}
-		getline(in, line); //$EndMeshFormat
+		// Node data
 		getline(in, line); //$Nodes Keywords
-
-		size_t n_nodes(0);
-		size_t n_elements(0);
-		while (line.find("$EndElements") == std::string::npos)
+		if (line.find("$Nodes") != std::string::npos)
 		{
-			// Node data
+			size_t n_nodes(0);
 			long id;
 			double x, y, z;
 			in >> n_nodes >> std::ws;
 			nodes.resize(n_nodes);
-			std::map<unsigned, unsigned> id_map;
 			for (size_t i = 0; i < n_nodes; i++) {
 				in >> id >> x >> y >> z >> std::ws;
 				id_map.insert(std::map<unsigned, unsigned>::value_type(id, i));
 				nodes[i] = new MeshLib::Node(x,y,z,id);
 			}
 			getline(in, line); // End Node keyword $EndNodes
+		}
 
-			// Element data
-			getline(in, line); // Element keyword $Elements
+		// Element data
+		if (line.find("$Elements") != std::string::npos)
+		{
+			size_t n_elements(0);
 			in >> n_elements >> std::ws; // number-of-elements
 			elements.reserve(n_elements);
-			unsigned idx, type, n_tags, dummy, mat_id;
 			for (size_t i = 0; i < n_elements; i++)
 			{
-				MeshLib::Element* elem (NULL);
-				std::vector<unsigned> node_ids;
-				std::vector<MeshLib::Node*> elem_nodes;
-				in >> idx >> type >> n_tags >> dummy >> mat_id;
-
-				// skip tags
-				for (size_t j = 2; j < n_tags; j++)
-					in >> dummy;
-
-				switch (type)
-				{
-				case 1: {
-					readNodeIDs(in, 2, node_ids, id_map);
-					// edge_nodes array will be deleted from Line object
-					MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
-					edge_nodes[0] = nodes[node_ids[0]];
-					edge_nodes[1] = nodes[node_ids[1]];
-					elem = new MeshLib::Line(edge_nodes, 0);
-					break;
-				}
-				case 2: {
-					readNodeIDs(in, 3, node_ids, id_map);
-					MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
-					tri_nodes[0] = nodes[node_ids[2]];
-					tri_nodes[1] = nodes[node_ids[1]];
-					tri_nodes[2] = nodes[node_ids[0]];
-					elem = new MeshLib::Tri(tri_nodes, mat_id);
-					break;
-				}
-				case 3: {
-					readNodeIDs(in, 4, node_ids, id_map);
-					MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
-					for (unsigned k(0); k < 4; k++)
-						quad_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Quad(quad_nodes, mat_id);
-					break;
-				}
-				case 4: {
-					readNodeIDs(in, 4, node_ids, id_map);
-					MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
-					for (unsigned k(0); k < 4; k++)
-						tet_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Tet(tet_nodes, mat_id);
-					break;
-				}
-				case 5: {
-					readNodeIDs(in, 8, node_ids, id_map);
-					MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
-					for (unsigned k(0); k < 8; k++)
-						hex_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Hex(hex_nodes, mat_id);
-					break;
-				}
-				case 6: {
-					readNodeIDs(in, 6, node_ids, id_map);
-					MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
-					for (unsigned k(0); k < 6; k++)
-						prism_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Prism(prism_nodes, mat_id);
-					break;
-				}
-				case 7: {
-					readNodeIDs(in, 5, node_ids, id_map);
-					MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
-					for (unsigned k(0); k < 5; k++)
-						pyramid_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Pyramid(pyramid_nodes, mat_id);
-					break;
-				}
-				case 15:
-					in >> dummy; // skip rest of line
-					continue;
-					break;
-				default:
-						WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
-				}
-				in >> std::ws;
+				MeshLib::Element* elem (readElement(in, nodes, id_map));
 
-				if (type > 0 && type < 8)
+				if (elem)
 					elements.push_back(elem);
 			}
+			getline(in, line); // END keyword
+		}
 
+		if (line.find("PhysicalNames") != std::string::npos)
+		{
+			size_t n_lines(0);
+			in >> n_lines >> std::ws; // number-of-lines
+			for (size_t i = 0; i < n_lines; i++)
+				getline(in, line);
 			getline(in, line); // END keyword
 		}
 	}
@@ -243,16 +189,91 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 void GMSHInterface::readNodeIDs(std::ifstream &in,
                                 unsigned n_nodes,
                                 std::vector<unsigned> &node_ids,
-                                std::map<unsigned, unsigned> &id_map)
+                                std::map<unsigned, unsigned> const& id_map)
 {
 	unsigned idx;
 	for (unsigned i = 0; i < n_nodes; i++)
 	{
 		in >> idx;
-		node_ids.push_back(id_map[idx]);
+		node_ids.push_back(id_map.at(idx));
 	}
 }
 
+MeshLib::Element* GMSHInterface::readElement(std::ifstream &in, std::vector<MeshLib::Node*> const& nodes, std::map<unsigned, unsigned> const& id_map)
+{	
+	unsigned idx, type, n_tags, dummy, mat_id;
+	MeshLib::Element* elem (nullptr);
+	std::vector<unsigned> node_ids;
+	std::vector<MeshLib::Node*> elem_nodes;
+	in >> idx >> type >> n_tags >> dummy >> mat_id;
+
+	// skip tags
+	for (size_t j = 2; j < n_tags; j++)
+		in >> dummy;
+
+	switch (type)
+	{
+	case 1: {
+		readNodeIDs(in, 2, node_ids, id_map);
+		// edge_nodes array will be deleted from Line object
+		MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
+		edge_nodes[0] = nodes[node_ids[0]];
+		edge_nodes[1] = nodes[node_ids[1]];
+		return new MeshLib::Line(edge_nodes, 0);
+	}
+	case 2: {
+		readNodeIDs(in, 3, node_ids, id_map);
+		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
+		tri_nodes[0] = nodes[node_ids[2]];
+		tri_nodes[1] = nodes[node_ids[1]];
+		tri_nodes[2] = nodes[node_ids[0]];
+		return new MeshLib::Tri(tri_nodes, mat_id);
+	}
+	case 3: {
+		readNodeIDs(in, 4, node_ids, id_map);
+		MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
+		for (unsigned k(0); k < 4; k++)
+			quad_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Quad(quad_nodes, mat_id);
+	}
+	case 4: {
+		readNodeIDs(in, 4, node_ids, id_map);
+		MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
+		for (unsigned k(0); k < 4; k++)
+			tet_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Tet(tet_nodes, mat_id);
+	}
+	case 5: {
+		readNodeIDs(in, 8, node_ids, id_map);
+		MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
+		for (unsigned k(0); k < 8; k++)
+			hex_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Hex(hex_nodes, mat_id);
+	}
+	case 6: {
+		readNodeIDs(in, 6, node_ids, id_map);
+		MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
+		for (unsigned k(0); k < 6; k++)
+			prism_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Prism(prism_nodes, mat_id);
+	}
+	case 7: {
+		readNodeIDs(in, 5, node_ids, id_map);
+		MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
+		for (unsigned k(0); k < 5; k++)
+			pyramid_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Pyramid(pyramid_nodes, mat_id);
+	}
+	case 15:
+		in >> dummy; // skip rest of line
+		break;
+	default:
+		WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
+	}
+	return nullptr;
+}
+
+
 bool GMSHInterface::write()
 {
 	_out << "// GMSH input file created by OpenGeoSys " << BaseLib::BuildInfo::ogs_version_and_persons;

FileIO/GMSHInterface.h

@@ -30,6 +30,8 @@
 
 namespace MeshLib {
 	class Mesh;
+	class Element;
+	class Node;
 }
 
 namespace FileIO
@@ -87,6 +89,9 @@ protected:
 	bool write();
 
 private:
+	/// Reads a mesh element from the input stream
+	static MeshLib::Element* readElement(std::ifstream &in, std::vector<MeshLib::Node*> const& nodes, std::map<unsigned, unsigned> const& id_map);
+
 	/**
 	 * 1. get and merge data from _geo_objs
 	 * 2. compute topological hierarchy
@@ -94,7 +99,7 @@ private:
 	 */
 	void writeGMSHInputFile(std::ostream & out);
 
-	static void readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector<unsigned> &node_ids, std::map<unsigned, unsigned> &id_map);
+	static void readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector<unsigned> &node_ids, std::map<unsigned, unsigned> const& id_map);
 
 	void writePoints(std::ostream& out) const;
 

FileIO/Legacy/OGSIOVer4.cpp

@@ -565,6 +565,39 @@ bool readGLIFileV4(const std::string& fname,
 		return false;
 }
 
+std::size_t writeTINSurfaces(std::ofstream &os, GeoLib::SurfaceVec const* sfcs_vec, std::size_t sfc_count, std::string const& path)
+{
+	const std::vector<GeoLib::Surface*>* sfcs (sfcs_vec->getVector());
+	for (std::size_t k(0); k < sfcs->size(); k++)
+	{
+		os << "#SURFACE" << "\n";
+		std::string sfc_name;
+		if (sfcs_vec->getNameOfElementByID (sfc_count, sfc_name)) {
+			os << "\t$NAME " << "\n" << "\t\t" << sfc_name << "\n";
+		} else {
+			os << "\t$NAME " << "\n" << "\t\t" << sfc_count << "\n";
+			sfc_name = std::to_string (sfc_count);
+		}
+		sfc_name += ".tin";
+		os << "\t$TIN" << "\n";
+		os << "\t\t" << sfc_name << "\n";
+		// create tin file
+		sfc_name = path + sfc_name;
+		std::ofstream tin_os (sfc_name.c_str());
+		GeoLib::Surface const& sfc (*(*sfcs)[k]);
+		const std::size_t n_tris (sfc.getNTriangles());
+		for (std::size_t l(0); l < n_tris; l++) {
+			GeoLib::Triangle const& tri (*(sfc[l]));
+		tin_os << l << " " << *(tri.getPoint(0)) << " " <<
+		*(tri.getPoint(1)) << " " << *(tri.getPoint(2)) <<
+		"\n";
+		}
+		tin_os.close();
+		sfc_count++;
+	}
+	return sfc_count;
+}
+
 void writeGLIFileV4 (const std::string& fname,
                      const std::string& geo_name,
                      const GeoLib::GEOObjects& geo)
@@ -603,19 +636,13 @@ void writeGLIFileV4 (const std::string& fname,
 			for (std::size_t j(0); j < (*plys)[k]->getNumberOfPoints(); j++)
 				os << "  " << ((*plys)[k])->getPointID(j) << "\n";
 		}
-
-		INFO("GeoLib::writeGLIFileV4(): write closed polylines as surfaces to file %s.",
-		     fname.c_str());
-		for (std::size_t k(0); k < plys->size(); k++)
-			if ((*plys)[k]->isClosed())
-			{
-				os << "#SURFACE" << "\n";
-				os << " $NAME " << "\n" << "  " << k << "\n"; //plys_vec->getNameOfElement ((*plys)[k]) << "\n";
-				os << " $TYPE " << "\n" << "  0" << "\n";
-				os << " $POLYLINES" << "\n" << "  " << k << "\n"; //plys_vec->getNameOfElement ((*plys)[k]) << "\n";
-			}
 	}
 
+	// writing surfaces as TIN files
+	const GeoLib::SurfaceVec* sfcs_vec (geo.getSurfaceVecObj (geo_name));
+	if (sfcs_vec) 
+		writeTINSurfaces(os, sfcs_vec, 0, BaseLib::extractPath(fname));
+
 	os << "#STOP" << "\n";
 	os.close ();
 }
@@ -706,37 +733,8 @@ void writeAllDataToGLIFileV4 (const std::string& fname, const GeoLib::GEOObjects
 	for (std::size_t j(0); j < geo_names.size(); j++)
 	{
 		const GeoLib::SurfaceVec* sfcs_vec (geo.getSurfaceVecObj (geo_names[j]));
-			if (sfcs_vec) {
-				const std::vector<GeoLib::Surface*>* sfcs (sfcs_vec->getVector());
-				for (std::size_t k(0); k < sfcs->size(); k++)
-				{
-					os << "#SURFACE" << "\n";
-					std::string sfc_name;
-					if (sfcs_vec->getNameOfElementByID (sfcs_cnt, sfc_name)) {
-						os << "\t$NAME " << "\n" << "\t\t" << sfc_name << "\n";
-					} else {
-						os << "\t$NAME " << "\n" << "\t\t" << sfcs_cnt << "\n";
-					sfc_name = std::to_string (sfcs_cnt);
-				}
-				sfc_name += ".tin";
-				os << "\t$TIN" << "\n";
-				os << "\t\t" << sfc_name << "\n";
-				// create tin file
-				sfc_name = path + sfc_name;
-				std::ofstream tin_os (sfc_name.c_str());
-				GeoLib::Surface const& sfc (*(*sfcs)[k]);
-				const std::size_t n_tris (sfc.getNTriangles());
-					for (std::size_t l(0); l < n_tris; l++) {
-						GeoLib::Triangle const& tri (*(sfc[l]));
-					tin_os << l << " " << *(tri.getPoint(0)) << " " <<
-					*(tri.getPoint(1)) << " " << *(tri.getPoint(2)) <<
-					"\n";
-				}
-				tin_os.close();
-
-				sfcs_cnt++;
-			}
-		}
+		if (sfcs_vec)
+			sfcs_cnt += writeTINSurfaces(os, sfcs_vec, sfcs_cnt, path);
 	}
 
 	os << "#STOP" << "\n";

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);

GeoLib/MinimalBoundingSphere.h

@@ -31,7 +31,7 @@ class MinimalBoundingSphere
 {
 public:
     /// Copy constructor
-    MinimalBoundingSphere(MinimalBoundingSphere const&) = default;
+    MinimalBoundingSphere(MinimalBoundingSphere const& sphere) = default;
     /// Point-Sphere
     MinimalBoundingSphere(GeoLib::Point const& p, double radius = std::numeric_limits<double>::epsilon());
     /// Bounding sphere using two points

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

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);
 

GeoLib/TemplateVec.h

@@ -82,6 +82,12 @@ public:
 	 */
 	std::string getName () const { return _name; }
 
+	/// Returns the begin of the name id mapping structure
+	NameIdMap::const_iterator getNameIDMapBegin() const { return _name_id_map->cbegin(); }
+	
+	/// Returns the end of the name id mapping structure
+	NameIdMap::const_iterator getNameIDMapEnd() const { return _name_id_map->cend(); }
+
 	/**
 	 * @return the number of data elements
 	 */

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;
 }
 

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;

Gui/DataView/GEOModels.cpp

@@ -46,14 +46,14 @@ void GEOModels::updateGeometry(const std::string &geo_name)
 	{
 		emit geoDataRemoved(_geoModel, geo_name, GeoLib::GEOTYPE::POINT);
 		this->_geoModel->removeGeoList(geo_name, GeoLib::GEOTYPE::POINT);
-		_geoModel->addPointList(QString::fromStdString(geo_name), points);
+		_geoModel->addPointList(QString::fromStdString(geo_name), *points);
 		emit geoDataAdded(_geoModel, geo_name, GeoLib::GEOTYPE::POINT);
 
 		if (lines)
 		{
 			emit geoDataRemoved(_geoModel, geo_name, GeoLib::GEOTYPE::POLYLINE);
 			this->_geoModel->removeGeoList(geo_name, GeoLib::GEOTYPE::POLYLINE);
-			_geoModel->addPolylineList(QString::fromStdString(geo_name), lines);
+			_geoModel->addPolylineList(QString::fromStdString(geo_name), *lines);
 			emit geoDataAdded(_geoModel, geo_name, GeoLib::GEOTYPE::POLYLINE);
 		}
 
@@ -61,7 +61,7 @@ void GEOModels::updateGeometry(const std::string &geo_name)
 		{
 			emit geoDataRemoved(_geoModel, geo_name, GeoLib::GEOTYPE::SURFACE);
 			this->_geoModel->removeGeoList(geo_name, GeoLib::GEOTYPE::SURFACE);
-			_geoModel->addSurfaceList(QString::fromStdString(geo_name), surfaces);
+			_geoModel->addSurfaceList(QString::fromStdString(geo_name), *surfaces);
 			emit geoDataAdded(_geoModel, geo_name, GeoLib::GEOTYPE::SURFACE);
 		}
 	}
@@ -85,7 +85,7 @@ void GEOModels::addPointVec( std::vector<GeoLib::Point*>* points,
                              double eps)
 {
 	GEOObjects::addPointVec(points, name, name_pnt_id_map, eps);
-	_geoModel->addPointList(QString::fromStdString(name), this->getPointVecObj(name));
+	_geoModel->addPointList(QString::fromStdString(name), *this->getPointVecObj(name));
 	emit geoDataAdded(_geoModel, name, GeoLib::GEOTYPE::POINT);
 }
 
@@ -141,7 +141,7 @@ void GEOModels::addPolylineVec( std::vector<GeoLib::Polyline*>* lines,
 	if (lines->empty())
 		return;
 
-	_geoModel->addPolylineList(QString::fromStdString(name), this->getPolylineVecObj(name));
+	_geoModel->addPolylineList(QString::fromStdString(name), *this->getPolylineVecObj(name));
 	emit geoDataAdded(_geoModel, name, GeoLib::GEOTYPE::POLYLINE);
 }
 
@@ -150,7 +150,7 @@ bool GEOModels::appendPolylineVec(const std::vector<GeoLib::Polyline*> &polyline
 {
 	bool ret (GeoLib::GEOObjects::appendPolylineVec (polylines, name));
 
-	this->_geoModel->appendPolylines(name, this->getPolylineVecObj(name));
+	this->_geoModel->appendPolylines(name, *this->getPolylineVecObj(name));
 	return ret;
 }
 
@@ -169,7 +169,7 @@ void GEOModels::addSurfaceVec( std::vector<GeoLib::Surface*>* surfaces,
 	if (surfaces->empty())
 		return;
 
-	_geoModel->addSurfaceList(QString::fromStdString(name), this->getSurfaceVecObj(name));
+	_geoModel->addSurfaceList(QString::fromStdString(name), *this->getSurfaceVecObj(name));
 	emit geoDataAdded(_geoModel, name, GeoLib::GEOTYPE::SURFACE);
 }
 
@@ -179,7 +179,7 @@ bool GEOModels::appendSurfaceVec(const std::vector<GeoLib::Surface*> &surfaces,
 	bool ret (GeoLib::GEOObjects::appendSurfaceVec (surfaces, name));
 
 	if (ret)
-		this->_geoModel->appendSurfaces(name, this->getSurfaceVecObj(name));
+		this->_geoModel->appendSurfaces(name, *this->getSurfaceVecObj(name));
 	else
 	{
 		std::vector<GeoLib::Surface*>* sfc = new std::vector<GeoLib::Surface*>;

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

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;

Gui/DataView/GeoTreeModel.cpp

@@ -36,9 +36,9 @@ GeoTreeModel::~GeoTreeModel()
 {
 }
 
-void GeoTreeModel::addPointList(QString geoName, const GeoLib::PointVec* pointVec)
+void GeoTreeModel::addPointList(QString geoName, GeoLib::PointVec const& pointVec)
 {
-	const std::vector<GeoLib::Point*>* points = pointVec->getVector();
+	const std::vector<GeoLib::Point*>* points = pointVec.getVector();
 
 	QList<QVariant> geoData;
 	geoData << QVariant(geoName) << "" << "" << "" << "";
@@ -56,26 +56,27 @@ void GeoTreeModel::addPointList(QString geoName, const GeoLib::PointVec* pointVe
 	for (size_t j = 0; j < nPoints; j++)
 	{
 		const GeoLib::Point &pnt(*(*points)[j]);
-		std::string pnt_name("");
-		pointVec->getNameOfElementByID(j, pnt_name);
 		QList<QVariant> pnt_data;
 		pnt_data.reserve(5);
 		pnt_data << static_cast<unsigned>(j)
 		         << QString::number(pnt[0], 'f')
 		         << QString::number(pnt[1], 'f')
 		         << QString::number(pnt[2], 'f')
-		         << QString::fromStdString(pnt_name);
+		         << "";
 		pointList->appendChild(new GeoTreeItem(pnt_data,
 		                                       pointList,
 		                                       static_cast<const GeoLib::Point*>(&pnt)));
 	}
 
+	for (auto pnt = pointVec.getNameIDMapBegin(); pnt != pointVec.getNameIDMapEnd(); ++pnt)
+		QVariant pnt_data (pointList->child(pnt->second)->setData(4, QString::fromStdString(pnt->first)));
+
 	INFO("Geometry \"%s\" built. %d points added.", geoName.toStdString().c_str(), nPoints);
 
 	reset();
 }
 
-void GeoTreeModel::addPolylineList(QString geoName, const GeoLib::PolylineVec* polylineVec)
+void GeoTreeModel::addPolylineList(QString geoName, GeoLib::PolylineVec const& polylineVec)
 {
 	int nLists = _rootItem->childCount();
 	TreeItem* geo(NULL);
@@ -91,7 +92,7 @@ void GeoTreeModel::addPolylineList(QString geoName, const GeoLib::PolylineVec* p
 		return;
 	}
 
-	const std::vector<GeoLib::Polyline*>* lines = polylineVec->getVector();
+	const std::vector<GeoLib::Polyline*>* lines = polylineVec.getVector();
 
 	QList<QVariant> plyData;
 	plyData << "Polylines" << "" << "" << "";
@@ -101,7 +102,7 @@ void GeoTreeModel::addPolylineList(QString geoName, const GeoLib::PolylineVec* p
 	reset();
 }
 
-void GeoTreeModel::appendPolylines(const std::string &name, const GeoLib::PolylineVec* polylineVec)
+void GeoTreeModel::appendPolylines(const std::string &name, GeoLib::PolylineVec const& polylineVec)
 {
 	for (size_t i = 0; i < _lists.size(); i++)
 	{
@@ -114,7 +115,7 @@ void GeoTreeModel::appendPolylines(const std::string &name, const GeoLib::Polyli
 				{
 					this->addChildren(parent, polylineVec,
 					                  parent->childCount(),
-					                  polylineVec->getVector()->size());
+					                  polylineVec.getVector()->size());
 					reset();
 					parent->vtkSource()->Modified();
 					return;
@@ -125,20 +126,17 @@ void GeoTreeModel::appendPolylines(const std::string &name, const GeoLib::Polyli
 }
 
 void GeoTreeModel::addChildren(GeoObjectListItem* plyList,
-                               const GeoLib::PolylineVec* polyline_vec,
+                               GeoLib::PolylineVec const& polyline_vec,
                                size_t start_index,
                                size_t end_index)
 {
-	const std::vector<GeoLib::Polyline*> lines = *(polyline_vec->getVector());
+	const std::vector<GeoLib::Polyline*> lines = *(polyline_vec.getVector());
 
 	for (size_t i = start_index; i < end_index; i++)
 	{
 		QList<QVariant> line_data;
 		line_data.reserve(4);
-		std::string ply_name("");
-		if (polyline_vec->getNameOfElementByID(i, ply_name))
-			line_data << "Line " + QString::number(i) << QString::fromStdString(ply_name) << "" << "";
-		else line_data << "Line " + QString::number(i) << "" << "" << "";
+		line_data << "Line " + QString::number(i) << "" << "" << "";
 
 		const GeoLib::Polyline &line(*(lines[i]));
 		GeoTreeItem* lineItem(new GeoTreeItem(line_data, plyList, &line));
@@ -158,10 +156,14 @@ void GeoTreeModel::addChildren(GeoObjectListItem* plyList,
 			lineItem->appendChild(new TreeItem(pnt_data, lineItem));
 		}
 	}
+
+	for (auto pnt = polyline_vec.getNameIDMapBegin(); pnt != polyline_vec.getNameIDMapEnd(); ++pnt)
+		QVariant pnt_data (plyList->child(pnt->second)->setData(4, QString::fromStdString(pnt->first)));
+
 	INFO("%d polylines added.", end_index - start_index);
 }
 
-void GeoTreeModel::addSurfaceList(QString geoName, const GeoLib::SurfaceVec* surfaceVec)
+void GeoTreeModel::addSurfaceList(QString geoName, GeoLib::SurfaceVec const& surfaceVec)
 {
 	int nLists = _rootItem->childCount();
 	TreeItem* geo(NULL);
@@ -177,7 +179,7 @@ void GeoTreeModel::addSurfaceList(QString geoName, const GeoLib::SurfaceVec* sur
 		return;
 	}
 
-	const std::vector<GeoLib::Surface*>* surfaces = surfaceVec->getVector();
+	const std::vector<GeoLib::Surface*>* surfaces = surfaceVec.getVector();
 
 	QList<QVariant> sfcData;
 	sfcData << "Surfaces" << "" << "" << "";
@@ -188,7 +190,7 @@ void GeoTreeModel::addSurfaceList(QString geoName, const GeoLib::SurfaceVec* sur
 	reset();
 }
 
-void GeoTreeModel::appendSurfaces(const std::string &name, GeoLib::SurfaceVec* surfaceVec)
+void GeoTreeModel::appendSurfaces(const std::string &name, GeoLib::SurfaceVec const& surfaceVec)
 {
 	for (size_t i = 0; i < _lists.size(); i++)
 	{
@@ -203,7 +205,7 @@ void GeoTreeModel::appendSurfaces(const std::string &name, GeoLib::SurfaceVec* s
 				{
 					this->addChildren(parent, surfaceVec,
 					                  parent->childCount(),
-					                  surfaceVec->getVector()->size());
+					                  surfaceVec.getVector()->size());
 					parent->vtkSource()->Modified();
 					reset();
 					return;
@@ -215,11 +217,11 @@ void GeoTreeModel::appendSurfaces(const std::string &name, GeoLib::SurfaceVec* s
 }
 
 void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
-                               const GeoLib::SurfaceVec* surface_vec,
+                               GeoLib::SurfaceVec const& surface_vec,
                                size_t start_index,
                                size_t end_index)
 {
-	const std::vector<GeoLib::Surface*>* surfaces = surface_vec->getVector();
+	const std::vector<GeoLib::Surface*>* surfaces = surface_vec.getVector();
 
 	const std::vector<GeoLib::Point*> &nodesVec(*((*surfaces)[start_index]->getPointVec()));
 	for (size_t i = start_index; i < end_index; i++)
@@ -227,7 +229,7 @@ void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
 		QList<QVariant> surface;
 		surface.reserve(4);
 		std::string sfc_name("");
-		surface_vec->getNameOfElementByID(i, sfc_name);
+		surface_vec.getNameOfElementByID(i, sfc_name);
 		surface << "Surface " + QString::number(i) << QString::fromStdString(sfc_name) <<
 		"" << "";
 

Gui/DataView/GeoTreeModel.h

@@ -51,15 +51,15 @@ public:
 	 * \param geoName Name of the new subtree. If no name is given a default name is assigned.
 	 * \param pointVec The list of points to be added as children of that subtree (no geometry can be added with a set of points!)
 	 */
-	void addPointList(QString geoName, const GeoLib::PointVec* pointVec);
+	void addPointList(QString geoName, GeoLib::PointVec const& pointVec);
 	/// Adds a subtree "Polylines" to an existing geometry with the given name.
-	void addPolylineList(QString geoName, const GeoLib::PolylineVec* polylineVec);
+	void addPolylineList(QString geoName, GeoLib::PolylineVec const& polylineVec);
 	/// Appends polylines to the "Polyline"-subtree
-	void appendPolylines(const std::string &name, const GeoLib::PolylineVec* polylineVec);
+	void appendPolylines(const std::string &name, GeoLib::PolylineVec const& polylineVec);
 	/// Adds a subtree "Surfaces" to an existing geometry with the given name.
-	void addSurfaceList(QString geoName, const GeoLib::SurfaceVec* surfaceVec);
+	void addSurfaceList(QString geoName, GeoLib::SurfaceVec const& surfaceVec);
 	/// Appends surfaces to the "Surface"-subtree
-	void appendSurfaces(const std::string &name, GeoLib::SurfaceVec* surfaceVec);
+	void appendSurfaces(const std::string &name, GeoLib::SurfaceVec const& surfaceVec);
 	/// Returns a list of all existing geometries.
 	const std::vector<GeoTreeItem*> &getLists() { return _lists; }
 	/**
@@ -84,13 +84,13 @@ public:
 private:
 	/// Adds children to the "Polylines" node
 	void addChildren(GeoObjectListItem* plyList,
-	                 const GeoLib::PolylineVec* polyline_vec,
+	                 GeoLib::PolylineVec const& polyline_vec,
 	                 std::size_t start_index,
 	                 std::size_t end_index);
 
 	/// Adds children to the "Surfaces" node
 	void addChildren(GeoObjectListItem* sfcList,
-	                 const GeoLib::SurfaceVec* surface_vec,
+	                 GeoLib::SurfaceVec const& surface_vec,
 	                 std::size_t start_index,
 	                 std::size_t end_index);
 

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;