sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-02-22 00:06:12 -05:00
Code Issues Packages Projects Releases Wiki Activity

update MDAL to 1.1.0 (#55457)

Browse Source 
update MDAL to 1.1.0
This commit is contained in:
Jan Caha 2023-12-01 12:15:26 +01:00 committed by GitHub
parent f141ce250f
commit 744a1691b8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 601 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

491
external/mdal/frmts/mdal_mike21.cpp vendored Normal file
View File

@ -0,0 +1,491 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2023 Lutra Consulting Ltd.
*/
#include <stddef.h>
#include <iosfwd>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cassert>
#include <limits>
#include <algorithm>
#include <regex>
#include "mdal_mike21.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
#include "mdal_logger.hpp"
#define DRIVER_NAME "Mike21"
// function to split using regex, by default split on whitespace characters
std::vector<std::string> regex_split( const std::string &input, const std::regex &split_regex = std::regex{"\\s+"} )
{
std::sregex_token_iterator iter( input.begin(), input.end(), split_regex, -1 );
std::sregex_token_iterator end;
return {iter, end};
}
static bool parse_vertex_id_gaps( std::map<size_t, size_t> &vertexIDtoIndex, size_t vertexIndex, size_t vertexID )
{
if ( vertexIndex == vertexID )
return false;
std::map<size_t, size_t>::iterator search = vertexIDtoIndex.find( vertexID );
if ( search != vertexIDtoIndex.end() )
{
MDAL::Log::warning( Warn_ElementNotUnique, DRIVER_NAME, "could not find vertex" );
return true;
}
vertexIDtoIndex[vertexID] = vertexIndex;
return false;
}
static void persist_native_index( std::vector<double> &arr, size_t nativeID, size_t ourId, size_t maxOurId )
{
if ( !arr.empty() || ( nativeID != ourId + 1 ) )
{
// we have gaps in face indexing
if ( arr.empty() )
{
arr.resize( maxOurId );
for ( size_t i = 0; i < ourId; ++i )
arr[i] = static_cast<double>( i + 1 );
}
arr[ourId] = static_cast<double>( nativeID );
}
}
MDAL::MeshMike21::MeshMike21( size_t faceVerticesMaximumCount,
const std::string &uri,
const std::map<size_t, size_t> vertexIDtoIndex )
: MemoryMesh( DRIVER_NAME,
faceVerticesMaximumCount,
uri )
, mVertexIDtoIndex( vertexIDtoIndex )
{
}
MDAL::MeshMike21::~MeshMike21() = default;
size_t MDAL::MeshMike21::vertexIndex( size_t vertexID ) const
{
auto ni2i = mVertexIDtoIndex.find( vertexID );
if ( ni2i != mVertexIDtoIndex.end() )
{
return ni2i->second; // convert from ID to index
}
return vertexID;
}
size_t MDAL::MeshMike21::maximumVertexId() const
{
size_t maxIndex = verticesCount() - 1;
if ( mVertexIDtoIndex.empty() )
return maxIndex;
else
{
// std::map is sorted!
size_t maxID = mVertexIDtoIndex.rbegin()->first;
return std::max( maxIndex, maxID );
}
}
MDAL::DriverMike21::DriverMike21( ):
Driver( DRIVER_NAME,
"Mike21 Mesh File",
"*.mesh",
Capability::ReadMesh | Capability::SaveMesh
)
{
}
MDAL::DriverMike21 *MDAL::DriverMike21::create()
{
return new DriverMike21();
}
MDAL::DriverMike21::~DriverMike21() = default;
bool MDAL::DriverMike21::canReadHeader( const std::string &line )
{
bool header2012 = std::regex_match( line, mRegexHeader2012 );
bool header2011 = std::regex_match( line, mRegexHeader2011 );
return header2011 || header2012;
}
bool MDAL::DriverMike21::canReadMesh( const std::string &uri )
{
std::ifstream in = MDAL::openInputFile( uri );
std::string line;
if ( !MDAL::getHeaderLine( in, line ) || !canReadHeader( line ) || !MDAL::contains( filters(), MDAL::fileExtension( uri ) ) )
{
return false;
}
return true;
}
size_t MDAL::DriverMike21::getVertexCount( const std::string &line )
{
auto matchResults = std::smatch{};
if ( std::regex_search( line, matchResults, mRegexHeader2012 ) )
{
if ( matchResults.size() > 4 )
{
return std::stoi( matchResults[3].str() );
}
}
if ( std::regex_search( line, matchResults, mRegexHeader2011 ) )
{
if ( matchResults.size() > 2 )
{
return std::stoi( matchResults[1].str() );
}
}
return 0;
}
std::string MDAL::DriverMike21::getCrs( const std::string &line )
{
auto matchResults = std::smatch{};
if ( std::regex_search( line, matchResults, mRegexHeader2012 ) )
{
if ( matchResults.size() > 5 )
{
return matchResults[4].str();
}
}
if ( std::regex_search( line, matchResults, mRegexHeader2011 ) )
{
if ( matchResults.size() > 3 )
{
return matchResults[2].str();
}
}
return "";
}
std::unique_ptr<MDAL::Mesh> MDAL::DriverMike21::load( const std::string &meshFile, const std::string & )
{
mMeshFile = meshFile;
MDAL::Log::resetLastStatus();
std::ifstream in = MDAL::openInputFile( meshFile );
std::string line;
if ( !std::getline( in, line ) || !canReadHeader( line ) )
{
MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), meshFile + " could not be opened" );
return nullptr;
}
std::string crs = getCrs( line );
size_t vertexCount = getVertexCount( line );
size_t faceCount = 0;
size_t maxVerticesPerFace = 2;
size_t lineNumber = 1;
while ( std::getline( in, line ) )
{
if ( lineNumber == vertexCount + 1 )
{
auto matchResults = std::smatch{};
if ( std::regex_search( line, matchResults, mRegexElementHeader ) )
{
if ( matchResults.size() >= 4 )
{
faceCount = MDAL::toSizeT( matchResults[1].str() );
maxVerticesPerFace = MDAL::toSizeT( matchResults[2].str() );
size_t meshType = MDAL::toSizeT( matchResults[3].str() );
if ( !( meshType == 21 || meshType == 25 ) )
{
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "unknow mesh type." );
return nullptr;
}
}
else
{
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "element header not in valid format." );
return nullptr;
}
}
else
{
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "element header not in valid format." );
return nullptr;
}
}
lineNumber++;
}
// number of lines in file does not match number of vertices and faces specifed in first and element line
if ( lineNumber > 2 + vertexCount + faceCount )
{
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "Number of lines in file does not fit with number of vertexes and faces specified." );
return nullptr;
}
in.clear();
in.seekg( 0, std::ios::beg );
Vertices vertices( vertexCount );
Faces faces( faceCount );
std::map<size_t, size_t> vertexIDtoIndex;
std::vector<double> vertexType( vertexCount );
std::vector<double> nativeVertexIds;
std::vector<double> nativeFaceIds;
size_t lastVertexID = 0;
size_t faceIndex = 0;
size_t vertexIndex = 0;
std::vector<std::string> chunks;
lineNumber = 0;
while ( std::getline( in, line ) )
{
if ( 0 < lineNumber && lineNumber < vertexCount + 1 )
{
chunks = regex_split( MDAL::trim( line ) );
if ( chunks.size() != 5 )
{
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "vertex line in invalid format." );
return nullptr;
}
size_t nodeID = toSizeT( chunks[0] );
if ( nodeID != 0 )
{
// specification of Mike21 does not state if vertexIDs need to continuos, expect that they might be not
// in the same way as in 2DM
if ( ( lastVertexID != 0 ) && ( nodeID <= lastVertexID ) )
{
// the algorithm requires that the file has points orderer by index
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "nodes are not ordered by index" );
return nullptr;
}
lastVertexID = nodeID;
}
// in case we have gaps/reorders in native indexes, store it
persist_native_index( nativeVertexIds, nodeID, vertexIndex, vertexCount );
parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID - 1 );
assert( vertexIndex < vertexCount );
Vertex &vertex = vertices[vertexIndex];
vertex.x = toDouble( chunks[1] );
vertex.y = toDouble( chunks[2] );
vertex.z = toDouble( chunks[3] );
vertexType[vertexIndex] = MDAL::toInt( chunks[4] );
vertexIndex++;
}
if ( vertexCount + 1 < lineNumber )
{
chunks = regex_split( MDAL::trim( line ) );
assert( faceIndex < faceCount );
size_t faceVertexCount = chunks.size() - 1;
// if the face should have 4 vertexes last chunk has value 0
// it actually means that there are only 3 vertexes
if ( faceVertexCount == 4 && chunks.size() == 5 )
{
if ( MDAL::toSizeT( chunks[4] ) == 0 )
faceVertexCount = faceVertexCount - 1;
}
assert( ( faceVertexCount == 3 ) || ( faceVertexCount == 4 ) );
if ( maxVerticesPerFace < faceVertexCount )
maxVerticesPerFace = faceVertexCount;
Face &face = faces[faceIndex];
face.resize( faceVertexCount );
// in case we have gaps/reorders in native indexes, store it
size_t nativeID = MDAL::toSizeT( chunks[0] );
persist_native_index( nativeFaceIds, nativeID, faceIndex, faceCount );
for ( size_t i = 0; i < faceVertexCount; ++i )
face[i] = MDAL::toSizeT( chunks[i + 1] ) - 1; // Mike21 is numbered from 1
faceIndex++;
}
lineNumber++;
}
for ( std::vector<Face>::iterator it = faces.begin(); it != faces.end(); ++it )
{
Face &face = *it;
for ( Face::size_type nd = 0; nd < face.size(); ++nd )
{
size_t nodeID = face[nd];
std::map<size_t, size_t>::iterator ni2i = vertexIDtoIndex.find( nodeID );
if ( ni2i != vertexIDtoIndex.end() )
{
face[nd] = ni2i->second; // convert from ID to index
}
else if ( vertices.size() < nodeID )
{
MDAL::Log::warning( MDAL_Status::Warn_ElementWithInvalidNode, name(), "found invalid node" );
}
}
}
// create the mesh and set the required data
std::unique_ptr< MeshMike21 > mesh(
new MeshMike21(
maxVerticesPerFace,
mMeshFile,
vertexIDtoIndex
)
);
mesh->setFaces( std::move( faces ) );
mesh->setVertices( std::move( vertices ) );
// Add Vertex Type
MDAL::addVertexScalarDatasetGroup( mesh.get(), vertexType, "VertexType" );
// Add Bed Elevation
MDAL::addBedElevationDatasetGroup( mesh.get(), mesh->vertices() );
if ( !nativeFaceIds.empty() )
MDAL::addFaceScalarDatasetGroup( mesh.get(), nativeFaceIds, "NativeFaceIds" );
if ( !nativeVertexIds.empty() )
MDAL::addVertexScalarDatasetGroup( mesh.get(), nativeVertexIds, "NativeVertexIds" );
mesh->setSourceCrs( crs );
return std::unique_ptr<Mesh>( mesh.release() );
}
void MDAL::DriverMike21::save( const std::string &fileName, const std::string &, MDAL::Mesh *mesh )
{
MDAL::Log::resetLastStatus();
std::ofstream file = MDAL::openOutputFile( fileName, std::ofstream::out );
if ( !file.is_open() )
{
MDAL::Log::error( MDAL_Status::Err_FailToWriteToDisk, name(), "Could not open file " + fileName );
}
std::string line = std::to_string( mesh->verticesCount() ) + " " + mesh->crs();
file << line << std::endl;
std::vector<double> vertexTypes;
std::shared_ptr<MDAL::DatasetGroup> vertexTypeDG = mesh->group( "VertexType" );
if ( vertexTypeDG )
{
vertexTypes.resize( mesh->verticesCount() );
auto d = vertexTypeDG->datasets[0];
d->scalarData( 0, mesh->verticesCount(), vertexTypes.data() );
}
// write vertices
std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();
double vertex[3];
for ( size_t i = 0; i < mesh->verticesCount(); ++i )
{
vertexIterator->next( 1, vertex );
line = "";
line.append( std::to_string( i + 1 ) );
for ( size_t j = 0; j < 2; ++j )
{
line.append( " " );
line.append( MDAL::coordinateToString( vertex[j] ) );
}
line.append( " " );
line.append( MDAL::doubleToString( vertex[2] ) );
line.append( " " );
if ( vertexTypes.size() == mesh->verticesCount() )
{
line.append( MDAL::doubleToString( vertexTypes.at( i ) ) );
}
else
{
line.append( MDAL::doubleToString( 0 ) );
}
file << line << std::endl;
}
//write element header line
size_t elementType = 0;
if ( mesh->faceVerticesMaximumCount() == 3 )
{
elementType = 21;
}
else if ( mesh->faceVerticesMaximumCount() == 4 )
{
elementType = 25;
}
line = std::to_string( mesh->facesCount() );
line.append( " " );
line.append( std::to_string( mesh->faceVerticesMaximumCount() ) );
line.append( " " );
line.append( std::to_string( elementType ) );
file << line << std::endl;
// write faces
std::vector<int> vertexIndices( mesh->faceVerticesMaximumCount() );
std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
for ( size_t i = 0; i < mesh->facesCount(); ++i )
{
int faceOffsets[1];
faceIterator->next( 1, faceOffsets, 4, vertexIndices.data() );
if ( faceOffsets[0] > 2 && faceOffsets[0] < 5 )
{
line = "";
line.append( std::to_string( i + 1 ) );
for ( int j = 0; j < faceOffsets[0]; ++j )
{
line.append( " " );
line.append( std::to_string( vertexIndices[j] + 1 ) );
}
// if face has 3 vertexes but the mesh as whole is marked as having
// 4 vertex at maximum, the last element should 0 - indicating no vertex there
if ( faceOffsets[0] == 3 && mesh->faceVerticesMaximumCount() == 4 )
{
line.append( " " );
line.append( "0" );
}
}
file << line << std::endl;
}
file.close();
}
std::string MDAL::DriverMike21::saveMeshOnFileSuffix() const
{
return "mesh";
}

104
external/mdal/frmts/mdal_mike21.hpp vendored Normal file
View File

@ -0,0 +1,104 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2023 Lutra Consulting Ltd.
*/
#ifndef MDAL_MIKE21_HPP
#define MDAL_MIKE21_HPP
#include <string>
#include <memory>
#include <regex>
#include "mdal_data_model.hpp"
#include "mdal_memory_data_model.hpp"
#include "mdal.h"
#include "mdal_driver.hpp"
namespace MDAL
{
class MeshMike21: public MemoryMesh
{
public:
MeshMike21( size_t faceVerticesMaximumCount,
const std::string &uri,
const std::map<size_t, size_t> vertexIDtoIndex
);
~MeshMike21() override;
/**
* Mike21 may supports gaps in the vertex indexing,
* but we use continuos array of vertices in MDAL
* \param vertexID internal index/ID of the vertex that native format uses
* \returns index of the vertex in the continuous array of vertices we returned by readVertices().
* For invalid vertexID it is returned index that is out of vertices array bounds.
*/
virtual size_t vertexIndex( size_t vertexID ) const;
/**
* Returns maximum vertex ID.
* For meshes without gaps in vertex indexing, it is vertex count - 1
*/
virtual size_t maximumVertexId() const;
private:
/**
* Mike21 might supports "gaps" in the mesh indexing
* Store only the indices that have different index and ID
* https://github.com/lutraconsulting/MDAL/issues/51
*/
std::map<size_t, size_t> mVertexIDtoIndex;
};
/**
* Mike21 format specification
* Text file format representing mesh vertices and faces
* Vertices are stored as - id x y z type_of_vertex (0 for water, 1 for land and above 1 for all other boundaries)
* The format supports triangles and quads
* Faces are stored as id 1 2 3 [4]
*
* full specification here: https://www.xmswiki.com/wiki/SMS:MIKE_21_*.mesh
*
* The format has two special lines in the file:
* First line - that can have format either:
* int int int string - these are type of bathymetry data, units of bathymetry data, number of vertices, CRS
* int string - these are number of vertices, CRS
* Line after the vertices (Element header line):
* int int int - number of faces, max number of vertexes per face, mesh type code
* [21 for faces with triangles only, 25 for faces with also quadrangular elements]
*
* The driver creates at least two Vertex Scalar datasets on the Mesh. "Bed Elevation" for Z coordinate
* and "VertexType" storing the type of vertex (mentioned above). Besides "NativeVertexIds" and "NativeFaceIds"
* maybe created if the IDS in input data are not continuous.
*/
class DriverMike21: public Driver
{
public:
DriverMike21();
~DriverMike21() override;
DriverMike21 *create() override;
int faceVerticesMaximumCount() const override {return 4;}
bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, const std::string &meshName = "" ) override;
void save( const std::string &fileName, const std::string &, Mesh *mesh ) override;
std::string saveMeshOnFileSuffix() const override;
private:
std::string mMeshFile;
// regex for header line in form of - integer string
const std::regex mRegexHeader2011 = std::regex( "(\\d+)\\s+(.+)(\\s+)?" );
// regex for header line in form of - integer integer integer string
const std::regex mRegexHeader2012 = std::regex( "(\\d+)\\s+(\\d+)\\s+(\\d+)\\s+(.+)(\\s+)?" );
// regex for element header line - integer integer integer(code with two number)
const std::regex mRegexElementHeader = std::regex( "(\\d+)\\s+(\\d)\\s+(\\d{2})(\\s+)?" );
bool canReadHeader( const std::string &line );
size_t getVertexCount( const std::string &line );
std::string getCrs( const std::string &line );
};
} // namespace MDAL
#endif //MDAL_MIKE21_HPP

2
external/mdal/mdal.cpp vendored
View File

@ -21,7 +21,7 @@ static const char *EMPTY_STR = "";
const char *MDAL_Version()
{
return "1.0.3";
return "1.1.0";
}
MDAL_Status MDAL_LastStatus()

3
external/mdal/mdal_driver_manager.cpp vendored
View File

@ -11,6 +11,7 @@
#include "frmts/mdal_binary_dat.hpp"
#include "frmts/mdal_selafin.hpp"
#include "frmts/mdal_esri_tin.hpp"
#include "frmts/mdal_mike21.hpp"
#include "frmts/mdal_dynamic_driver.hpp"
#include "mdal_utils.hpp"
@ -262,6 +263,8 @@ MDAL::DriverManager::DriverManager()
mDrivers.push_back( std::make_shared<MDAL::DriverXdmf>() );
#endif
mDrivers.push_back( std::make_shared<MDAL::DriverMike21>() );
loadDynamicDrivers();
}

2
src/providers/mdal/CMakeLists.txt
View File

@ -48,6 +48,7 @@ if (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_esri_tin.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_xms_tin.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ply.cpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_mike21.cpp
${CMAKE_SOURCE_DIR}/external/mdal/3rdparty/libplyxx/libplyxx.cpp
)
@ -67,6 +68,7 @@ if (WITH_INTERNAL_MDAL)
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_esri_tin.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_xms_tin.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_ply.hpp
${CMAKE_SOURCE_DIR}/external/mdal/frmts/mdal_mike21.hpp
${CMAKE_SOURCE_DIR}/external/mdal/3rdparty/libplyxx.h
${CMAKE_SOURCE_DIR}/external/mdal/3rdparty/libplyxx/libplyxx_internal.h
${CMAKE_SOURCE_DIR}/external/mdal/3rdparty/textio.h