OGS
RasterToMesh.cpp
Go to the documentation of this file.
1
11#include "RasterToMesh.h"
12
13#include <vtkCell.h>
14#include <vtkCellData.h>
15#include <vtkDataArray.h>
16#include <vtkImageData.h>
17#include <vtkPointData.h>
18#include <vtkSmartPointer.h>
19
20#include "BaseLib/Logging.h"
21#include "MeshGenerator.h"
23#include "MeshLib/Mesh.h"
25#include "MeshLib/Node.h"
28
29namespace MeshToolsLib
30{
31std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert(
32 GeoLib::Raster const& raster,
33 MeshLib::MeshElemType elem_type,
34 MeshLib::UseIntensityAs intensity_type,
35 std::string const& array_name)
36{
37 return convert(raster.data(),
38 raster.getHeader(),
39 elem_type,
40 intensity_type,
41 array_name);
42}
43
44std::unique_ptr<MeshLib::Mesh> RasterToMesh::convert(
45 double const* const img,
46 GeoLib::RasterHeader const& header,
47 MeshLib::MeshElemType elem_type,
48 MeshLib::UseIntensityAs intensity_type,
49 std::string const& array_name)
50{
51 if ((elem_type != MeshLib::MeshElemType::TRIANGLE) &&
52 (elem_type != MeshLib::MeshElemType::QUAD) &&
53 (elem_type != MeshLib::MeshElemType::HEXAHEDRON) &&
54 (elem_type != MeshLib::MeshElemType::PRISM))
55 {
56 ERR("Invalid Mesh Element Type.");
57 return nullptr;
58 }
59
60 if (((elem_type == MeshLib::MeshElemType::TRIANGLE) ||
61 (elem_type == MeshLib::MeshElemType::QUAD)) &&
62 header.n_depth != 1)
63 {
64 ERR("Triangle or Quad elements cannot be used to construct meshes from "
65 "3D rasters.");
66 return nullptr;
67 }
68
69 if (intensity_type == MeshLib::UseIntensityAs::ELEVATION &&
70 ((elem_type == MeshLib::MeshElemType::PRISM) ||
72 {
73 ERR("Elevation mapping can only be performed for 2D meshes.");
74 return nullptr;
75 }
76
77 std::unique_ptr<MeshLib::Mesh> mesh;
78 if (elem_type == MeshLib::MeshElemType::TRIANGLE)
79 {
81 header.n_cols,
82 header.n_rows,
83 header.cell_size,
84 header.origin,
85 "RasterDataMesh"));
86 }
87 else if (elem_type == MeshLib::MeshElemType::QUAD)
88 {
90 header.n_cols,
91 header.n_rows,
92 header.cell_size,
93 header.origin,
94 "RasterDataMesh"));
95 }
96 else if (elem_type == MeshLib::MeshElemType::PRISM)
97 {
99 header.n_cols,
100 header.n_rows,
101 header.n_depth,
102 header.cell_size,
103 header.origin,
104 "RasterDataMesh"));
105 }
106 else if (elem_type == MeshLib::MeshElemType::HEXAHEDRON)
107 {
109 header.n_cols,
110 header.n_rows,
111 header.n_depth,
112 header.cell_size,
113 header.origin,
114 "RasterDataMesh"));
115 }
116
117 std::unique_ptr<MeshLib::Mesh> new_mesh;
118 std::vector<std::size_t> elements_to_remove;
119 if (intensity_type == MeshLib::UseIntensityAs::ELEVATION)
120 {
121 std::vector<MeshLib::Node*> const& nodes(mesh->getNodes());
122 std::vector<MeshLib::Element*> const& elems(mesh->getElements());
123 std::size_t const n_nodes(elems[0]->getNumberOfNodes());
124 bool const double_idx =
125 (elem_type == MeshLib::MeshElemType::TRIANGLE) ||
126 (elem_type == MeshLib::MeshElemType::PRISM);
127 std::size_t const m = (double_idx) ? 2 : 1;
128 for (std::size_t k = 0; k < header.n_depth; k++)
129 {
130 std::size_t const layer_idx = (k * header.n_rows * header.n_cols);
131 for (std::size_t i = 0; i < header.n_cols; i++)
132 {
133 std::size_t const idx(i * header.n_rows + layer_idx);
134 for (std::size_t j = 0; j < header.n_rows; j++)
135 {
136 double const val(img[idx + j]);
137 if (val == header.no_data)
138 {
139 elements_to_remove.push_back(m * (idx + j));
140 if (double_idx)
141 {
142 elements_to_remove.push_back(m * (idx + j) + 1);
143 }
144 continue;
145 }
146 for (std::size_t n = 0; n < n_nodes; ++n)
147 {
148 (*(nodes[getNodeIndex(*elems[m * (idx + j)], n)]))[2] =
149 val;
150 if (double_idx)
151 {
152 (*(nodes[getNodeIndex(*elems[m * (idx + j) + 1],
153 n)]))[2] = val;
154 }
155 }
156 }
157 }
158 }
159 }
160 else
161 {
162 MeshLib::Properties& properties = mesh->getProperties();
163 MeshLib::ElementSearch ex(*mesh);
164 if (array_name == "MaterialIDs")
165 {
166 auto* const prop_vec = properties.createNewPropertyVector<int>(
167 array_name, MeshLib::MeshItemType::Cell, 1);
168 fillPropertyVector<int>(*prop_vec, img, header, elem_type);
169 ex.searchByPropertyValue<int>(array_name,
170 static_cast<int>(header.no_data));
171 }
172 else
173 {
174 auto* const prop_vec = properties.createNewPropertyVector<double>(
175 array_name, MeshLib::MeshItemType::Cell, 1);
176 fillPropertyVector<double>(*prop_vec, img, header, elem_type);
177 ex.searchByPropertyValue<double>(array_name, header.no_data);
178 }
179 elements_to_remove = ex.getSearchedElementIDs();
180 }
181 if (!elements_to_remove.empty())
182 {
183 new_mesh.reset(MeshToolsLib::removeElements(
184 *mesh, elements_to_remove, mesh->getName()));
185 }
186 else
187 {
188 new_mesh = std::move(mesh);
189 }
190
191 if (intensity_type == MeshLib::UseIntensityAs::NONE)
192 {
193 new_mesh->getProperties().removePropertyVector(array_name);
194 }
195
196 return new_mesh;
197}
198
199} // namespace MeshToolsLib
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the Mesh class.
Definition of the Node class.
Class Raster is used for managing raster data.
Definition Raster.h:49
double const * data() const
Definition Raster.h:111
RasterHeader const & getHeader() const
Returns the complete header information.
Definition Raster.h:85
Element search class.
const std::vector< std::size_t > & getSearchedElementIDs() const
return marked elements
std::size_t searchByPropertyValue(std::string const &property_name, PROPERTY_TYPE const property_value)
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:36
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
static std::unique_ptr< MeshLib::Mesh > convert(GeoLib::Raster const &raster, MeshLib::MeshElemType elem_type, MeshLib::UseIntensityAs intensity_type, std::string const &array_name="Colour")
UseIntensityAs
Selection of possible interpretations for intensities.
Definition MeshEnums.h:83
MeshElemType
Types of mesh elements supported by OpenGeoSys. Values are from VTKCellType enum.
Definition MeshEnums.h:27
MeshLib::Mesh * generateRegularQuadMesh(const BaseLib::ISubdivision &div_x, const BaseLib::ISubdivision &div_y, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")
MeshLib::Mesh * generateRegularPrismMesh(const double x_length, const double y_length, const double z_length, const std::size_t x_subdivision, const std::size_t y_subdivision, const std::size_t z_subdivision, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")
MeshLib::Mesh * generateRegularTriMesh(const BaseLib::ISubdivision &div_x, const BaseLib::ISubdivision &div_y, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")
MeshLib::Mesh * generateRegularHexMesh(const BaseLib::ISubdivision &div_x, const BaseLib::ISubdivision &div_y, const BaseLib::ISubdivision &div_z, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)
Contains the relevant information when storing a geoscientific raster data.
Definition Raster.h:28
std::size_t n_depth
Definition Raster.h:31
MathLib::Point3d origin
Definition Raster.h:32
std::size_t n_cols
Definition Raster.h:29
std::size_t n_rows
Definition Raster.h:30