OGS
MeshMapping.cpp File Reference

Detailed Description

2012/03/07 KR Initial implementation

Definition in file MeshMapping.cpp.

#include <tclap/CmdLine.h>
#include <fstream>
#include <memory>
#include <string>
#include "BaseLib/FileTools.h"
#include "BaseLib/MPI.h"
#include "GeoLib/AABB.h"
#include "GeoLib/IO/AsciiRasterInterface.h"
#include "GeoLib/Raster.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/MathTools.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/IO/writeMeshToFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/MeshElementGrid.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshEditing/ProjectPointOnMesh.h"
#include "MeshToolsLib/MeshGenerators/MeshLayerMapper.h"
Include dependency graph for MeshMapping.cpp:

Go to the source code of this file.

Functions

double getClosestPointElevation (MeshLib::Node const &p, std::vector< MeshLib::Node * > const &nodes, double const &max_dist)
 
int main (int argc, char *argv[])
 

Function Documentation

◆ getClosestPointElevation()

double getClosestPointElevation ( MeshLib::Node const & p,
std::vector< MeshLib::Node * > const & nodes,
double const & max_dist )

Definition at line 33 of file MeshMapping.cpp.

36{
37 double sqr_shortest_dist(max_dist);
38 double elevation(p[2]);
39 for (MeshLib::Node* node : nodes)
40 {
41 double sqr_dist = (p[0] - (*node)[0]) * (p[0] - (*node)[0]) +
42 (p[1] - (*node)[1]) * (p[1] - (*node)[1]);
43 if (sqr_dist < sqr_shortest_dist)
44 {
45 sqr_shortest_dist = sqr_dist;
46 elevation = (*node)[2];
47 }
48 }
49 return elevation;
50}

Referenced by main().

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 52 of file MeshMapping.cpp.

53{
54 BaseLib::MPI::Setup mpi_setup(argc, argv);
55 TCLAP::CmdLine cmd(
56 "Changes the elevation of 2D mesh nodes based on either raster data or "
57 "another 2D mesh. In addition, a low pass filter can be applied to "
58 "node elevation based connected nodes.\n\n"
59 "OpenGeoSys-6 software, version " +
61 ".\n"
62 "Copyright (c) 2012-2024, OpenGeoSys Community "
63 "(http://www.opengeosys.org)",
65 TCLAP::SwitchArg lowpass_arg(
66 "", "lowpass",
67 "Applies a lowpass filter to elevation over connected nodes.", false);
68 cmd.add(lowpass_arg);
69 TCLAP::ValueArg<double> map_static_arg(
70 "s", "static",
71 "Static elevation to map the input file to. This can be combined with "
72 "mapping based on rasters or other meshes to deal with locations where "
73 "no corresponding data exists.",
74 false, 0, "number");
75 cmd.add(map_static_arg);
76 TCLAP::ValueArg<double> max_dist_arg(
77 "d", "distance",
78 "Maximum distance to search for mesh nodes if there is no "
79 "corresponding data for input mesh nodes on the mesh it should be "
80 "mapped on. (Default value: 1)",
81 false, 1, "number");
82 cmd.add(max_dist_arg);
83 TCLAP::ValueArg<std::string> map_mesh_arg(
84 "m", "mesh", "2D mesh file (*.vtu) to map the input file on.", false,
85 "", "string");
86 cmd.add(map_mesh_arg);
87 TCLAP::ValueArg<std::string> map_raster_arg(
88 "r", "raster",
89 "Raster file (*.asc *.grd *.xyz) to map the input file on.", false, "",
90 "string");
91 cmd.add(map_raster_arg);
92 TCLAP::ValueArg<std::string> output_arg(
93 "o", "output", "Output mesh file (*.vtu)", true, "", "string");
94 cmd.add(output_arg);
95 TCLAP::ValueArg<std::string> input_arg(
96 "i", "input", "Input mesh file (*.vtu, *.msh)", true, "", "string");
97 cmd.add(input_arg);
98 cmd.parse(argc, argv);
99
100 std::unique_ptr<MeshLib::Mesh> mesh(
101 MeshLib::IO::readMeshFromFile(input_arg.getValue()));
102 if (mesh == nullptr)
103 {
104 ERR("Error reading mesh file.");
105 return EXIT_FAILURE;
106 }
107
108 if (!(map_static_arg.isSet() || map_raster_arg.isSet() ||
109 map_mesh_arg.isSet()))
110 {
111 ERR("Nothing to do. Please choose mapping based on a raster or mesh "
112 "file, or to a static value.");
113 return EXIT_FAILURE;
114 }
115
116 if (map_raster_arg.isSet() && map_mesh_arg.isSet())
117 {
118 ERR("Please select mapping based on *either* a mesh or a raster file.");
119 return EXIT_FAILURE;
120 }
121
122 // Maps the elevation of mesh nodes to static value
123 if (map_static_arg.isSet())
124 {
126 *mesh, map_static_arg.getValue());
127 }
128
129 // Maps the elevation of mesh nodes according to raster
130 if (map_raster_arg.isSet())
131 {
132 std::string const raster_path = map_raster_arg.getValue();
133 std::ifstream file_stream(raster_path, std::ifstream::in);
134 if (!file_stream.good())
135 {
136 ERR("Opening raster file {} failed.", raster_path);
137 return EXIT_FAILURE;
138 }
139 file_stream.close();
140
141 std::unique_ptr<GeoLib::Raster> const raster(
143 MeshToolsLib::MeshLayerMapper::layerMapping(*mesh, *raster, 0, true);
144 }
145
146 // Maps the elevation of mesh nodes according to a ground truth mesh
147 if (map_mesh_arg.isSet())
148 {
149 std::unique_ptr<MeshLib::Mesh> ground_truth(
150 MeshLib::IO::readMeshFromFile(map_mesh_arg.getValue()));
151 if (ground_truth == nullptr)
152 {
153 ERR("Error reading mesh file.");
154 return EXIT_FAILURE;
155 }
156
157 std::vector<MeshLib::Node*> const& nodes = mesh->getNodes();
158 MeshLib::MeshElementGrid const grid(*ground_truth);
159 auto const edgeLengths = minMaxEdgeLength(mesh->getElements());
160 double const max_edge = edgeLengths.second;
161 double const max_dist(pow(max_dist_arg.getValue(), 2));
162
163 for (MeshLib::Node* node : nodes)
164 {
165 MathLib::Point3d min_vol{{(*node)[0] - max_edge,
166 (*node)[1] - max_edge,
167 -std::numeric_limits<double>::max()}};
168 MathLib::Point3d max_vol{{(*node)[0] + max_edge,
169 (*node)[1] + max_edge,
170 std::numeric_limits<double>::max()}};
171 std::vector<const MeshLib::Element*> const& elems =
172 grid.getElementsInVolume(min_vol, max_vol);
173 auto const* element =
175 *node);
176 if (element != nullptr)
177 {
179 *element, *node);
180 }
181 else
182 {
183 (*node)[2] = getClosestPointElevation(
184 *node, ground_truth->getNodes(), max_dist);
185 }
186 }
187 }
188
189 // a simple lowpass filter for the elevation of mesh nodes using the
190 // elevation of each node weighted by 2 and the elevation of each connected
191 // node weighted by 1
192 if (lowpass_arg.isSet())
193 {
194 INFO("lowpass");
195 const std::size_t nNodes(mesh->getNumberOfNodes());
196 std::vector<MeshLib::Node*> nodes(mesh->getNodes());
197
198 std::vector<double> elevation(nNodes);
199 for (std::size_t i = 0; i < nNodes; i++)
200 {
201 elevation[i] = (*nodes[i])[2];
202 }
203
204 auto const& connections =
206 for (std::size_t i = 0; i < nNodes; i++)
207 {
208 auto const& conn_nodes(connections[nodes[i]->getID()]);
209 const unsigned nConnNodes(conn_nodes.size());
210 elevation[i] = (2 * (*nodes[i])[2]);
211 for (std::size_t j = 0; j < nConnNodes; ++j)
212 {
213 elevation[i] += (*conn_nodes[j])[2];
214 }
215 elevation[i] /= (nConnNodes + 2);
216 }
217
218 for (std::size_t i = 0; i < nNodes; i++)
219 {
220 (*nodes[i])[2] = elevation[i];
221 }
222 }
223
224 if (MeshLib::IO::writeMeshToFile(*mesh, output_arg.getValue()) != 0)
225 {
226 return EXIT_FAILURE;
227 }
228
229 INFO("Result successfully written.");
230 return EXIT_SUCCESS;
231}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
double getClosestPointElevation(MeshLib::Node const &p, std::vector< MeshLib::Node * > const &nodes, double const &max_dist)
static GeoLib::Raster * readRaster(std::string const &fname)
static bool layerMapping(MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double nodata_replacement=0.0, bool const ignore_nodata=false)
static bool mapToStaticValue(MeshLib::Mesh const &mesh, double value)
Maps the elevation of all mesh nodes to the specified static value.
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, std::set< std::string > variable_output_names)
std::vector< std::vector< Node * > > calculateNodesConnectedByElements(Mesh const &mesh)
Definition Mesh.cpp:308
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.
Definition Mesh.cpp:189
double getElevation(MeshLib::Element const &element, MathLib::Point3d const &node)
MeshLib::Element const * getProjectedElement(std::vector< const MeshLib::Element * > const &elements, MathLib::Point3d const &node)

References MeshLib::calculateNodesConnectedByElements(), ERR(), getClosestPointElevation(), MeshLib::MeshElementGrid::getElementsInVolume(), MeshToolsLib::ProjectPointOnMesh::getElevation(), MeshToolsLib::ProjectPointOnMesh::getProjectedElement(), INFO(), MeshToolsLib::MeshLayerMapper::layerMapping(), MeshToolsLib::MeshLayerMapper::mapToStaticValue(), GitInfoLib::GitInfo::ogs_version, MeshLib::IO::readMeshFromFile(), FileIO::AsciiRasterInterface::readRaster(), and MeshLib::IO::writeMeshToFile().