QGIS API Documentation 3.41.0-Master (3440c17df1d)
Loading...
Searching...
No Matches
qgsmeshtriangulation.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsmeshtriangulation.cpp
3 -----------------
4 begin : August 9th, 2020
5 copyright : (C) 2020 by Vincent Cloarec
6 email : vcloarec at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19#include "moc_qgsmeshtriangulation.cpp"
21#include "qgscurve.h"
22#include "qgscurvepolygon.h"
23#include "qgsmultisurface.h"
24#include "qgsmulticurve.h"
25#include "qgsfeedback.h"
26#include "qgslogger.h"
27#include "qgsmesheditor.h"
28#include "qgsfeature.h"
29#include "qgsfeatureiterator.h"
30
32{
33 mTriangulation.reset( new QgsDualEdgeTriangulation() );
34}
35
36
38
39bool QgsMeshTriangulation::addVertices( QgsFeatureIterator &vertexFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback, long featureCount )
40{
41 if ( feedback )
42 feedback->setProgress( 0 );
43
44 QgsFeature feat;
45 long i = 0;
46 while ( vertexFeatureIterator.nextFeature( feat ) )
47 {
48 if ( feedback )
49 {
50 if ( feedback->isCanceled() )
51 break;
52
53 feedback->setProgress( 100 * i / featureCount );
54 i++;
55 }
56
57 addVerticesFromFeature( feat, valueAttribute, transform, feedback );
58 }
59
60 return true;
61}
62
63bool QgsMeshTriangulation::addBreakLines( QgsFeatureIterator &lineFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback, long featureCount )
64{
65 if ( feedback )
66 feedback->setProgress( 0 );
67
68 QgsFeature feat;
69 long i = 0;
70 while ( lineFeatureIterator.nextFeature( feat ) )
71 {
72 if ( feedback )
73 {
74 if ( feedback->isCanceled() )
75 break;
76
77 feedback->setProgress( 100 * i / featureCount );
78 i++;
79 }
80
81 Qgis::GeometryType geomType = feat.geometry().type();
82 switch ( geomType )
83 {
85 addVerticesFromFeature( feat, valueAttribute, transform, feedback );
86 break;
89 addBreakLinesFromFeature( feat, valueAttribute, transform, feedback );
90 break;
91 default:
92 break;
93 }
94 }
95
96 return true;
97}
98
100{
101 return mTriangulation->addPoint( vertex );
102}
103
105{
106 return mTriangulation->triangulationToMesh( feedback );
107}
108
113
114void QgsMeshTriangulation::addVerticesFromFeature( const QgsFeature &feature, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback )
115{
116 QgsGeometry geom = feature.geometry();
117 try
118 {
119 geom.transform( transform, Qgis::TransformDirection::Forward, true );
120 }
121 catch ( QgsCsException &cse )
122 {
123 Q_UNUSED( cse )
124 QgsDebugMsgLevel( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ), 4 );
125 return;
126 }
127
129
130 double value = 0;
131 if ( valueAttribute >= 0 )
132 value = feature.attribute( valueAttribute ).toDouble();
133
134 while ( vit != geom.vertices_end() )
135 {
136 if ( feedback && feedback->isCanceled() )
137 break;
138 if ( valueAttribute < 0 )
139 mTriangulation->addPoint( *vit );
140 else
141 {
142 mTriangulation->addPoint( QgsPoint( Qgis::WkbType::PointZ, ( *vit ).x(), ( *vit ).y(), value ) );
143 }
144 ++vit;
145 }
146}
147
148void QgsMeshTriangulation::addBreakLinesFromFeature( const QgsFeature &feature, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback )
149{
150 double valueOnVertex = 0;
151 if ( valueAttribute >= 0 )
152 valueOnVertex = feature.attribute( valueAttribute ).toDouble();
153
154 //from QgsTinInterpolator::insertData()
155 std::vector<const QgsCurve *> curves;
156 QgsGeometry geom = feature.geometry();
157 try
158 {
159 geom.transform( transform, Qgis::TransformDirection::Forward, true );
160 }
161 catch ( QgsCsException &cse )
162 {
163 Q_UNUSED( cse )
164 QgsDebugMsgLevel( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ), 4 );
165 return;
166 }
167
169 {
170 std::vector< const QgsCurvePolygon * > polygons;
171 if ( geom.isMultipart() )
172 {
173 const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( geom.constGet() );
174 for ( int i = 0; i < ms->numGeometries(); ++i )
175 {
176 polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( ms->geometryN( i ) ) );
177 }
178 }
179 else
180 {
181 polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( geom.constGet() ) );
182 }
183
184 for ( const QgsCurvePolygon *polygon : polygons )
185 {
186 if ( feedback && feedback->isCanceled() )
187 break;
188 if ( !polygon )
189 continue;
190
191 if ( polygon->exteriorRing() )
192 curves.emplace_back( polygon->exteriorRing() );
193
194 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
195 {
196 if ( feedback && feedback->isCanceled() )
197 break;
198 curves.emplace_back( polygon->interiorRing( i ) );
199 }
200 }
201 }
202 else
203 {
204 if ( geom.isMultipart() )
205 {
206 const QgsMultiCurve *mc = qgsgeometry_cast< const QgsMultiCurve * >( geom.constGet() );
207 for ( int i = 0; i < mc->numGeometries(); ++i )
208 {
209 if ( feedback && feedback->isCanceled() )
210 break;
211 curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( mc->geometryN( i ) ) );
212 }
213 }
214 else
215 {
216 curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( geom.constGet() ) );
217 }
218 }
219
220 for ( const QgsCurve *curve : curves )
221 {
222 if ( !curve )
223 continue;
224
225 if ( feedback && feedback->isCanceled() )
226 break;
227
228 QgsPointSequence linePoints;
229 curve->points( linePoints );
230 bool hasZ = curve->is3D();
231 if ( valueAttribute >= 0 )
232 for ( int i = 0; i < linePoints.count(); ++i )
233 {
234 if ( feedback && feedback->isCanceled() )
235 break;
236 if ( hasZ )
237 linePoints[i].setZ( valueOnVertex );
238 else
239 {
240 const QgsPoint &point = linePoints.at( i );
241 linePoints[i] = QgsPoint( point.x(), point.y(), valueOnVertex );
242 }
243 }
244
245 mTriangulation->addLine( linePoints, QgsInterpolator::SourceBreakLines );
246 }
247}
248
249QgsMeshZValueDatasetGroup::QgsMeshZValueDatasetGroup( const QString &datasetGroupName, const QgsMesh &mesh ):
250 QgsMeshDatasetGroup( datasetGroupName, QgsMeshDatasetGroupMetadata::DataOnVertices )
251{
252 mDataset = std::make_unique< QgsMeshZValueDataset >( mesh );
253}
254
259
261{
262 if ( datasetIndex != 0 )
263 return QgsMeshDatasetMetadata();
264
265 return mDataset->metadata();
266}
267
269
271{
272 if ( index != 0 )
273 return nullptr;
274
275 return mDataset.get();
276}
277
278QDomElement QgsMeshZValueDatasetGroup::writeXml( QDomDocument &doc, const QgsReadWriteContext &context ) const
279{
280 Q_UNUSED( doc );
281 Q_UNUSED( context );
282 return QDomElement();
283}
284
286{
287 for ( const QgsMeshVertex &vertex : mesh.vertices )
288 {
289 if ( vertex.z() < mZMinimum )
290 mZMinimum = vertex.z();
291 if ( vertex.z() > mZMaximum )
292 mZMaximum = vertex.z();
293 }
294}
295
297{
298 if ( valueIndex < 0 || valueIndex >= mMesh.vertexCount() )
299 return QgsMeshDatasetValue();
300
301 return QgsMeshDatasetValue( mMesh.vertex( valueIndex ).z() );
302}
303
304QgsMeshDataBlock QgsMeshZValueDataset::datasetValues( bool isScalar, int valueIndex, int count ) const
305{
306 Q_UNUSED( isScalar )
308 QVector<double> zValues( count );
309 for ( int i = valueIndex; i < valueIndex + count; ++i )
310 zValues[i - valueIndex] = mMesh.vertex( i ).z();
311 ret.setValues( zValues );
312 return ret;
313}
314
316{
317 Q_UNUSED( faceIndex );
318 Q_UNUSED( count );
320 ret.setValid( true );
321 return ret;
322}
323
324bool QgsMeshZValueDataset::isActive( int faceIndex ) const
325{
326 return ( faceIndex > 0 && faceIndex < mMesh.faceCount() );
327}
328
330{
331 return QgsMeshDatasetMetadata( 0, true, mZMinimum, mZMaximum, 0 );
332}
333
335{
336 return mMesh.vertexCount();
337}
338
340
342{
343 return QObject::tr( "Delaunay triangulation" );
344}
345
346QgsTopologicalMesh::Changes QgsMeshEditingDelaunayTriangulation::apply( QgsMeshEditor *meshEditor )
347{
348 //use only vertices that are on boundary or free, if boundary
349 QList<int> vertexIndextoTriangulate;
350
351 QList<int> removedVerticesFromTriangulation;
352
353 for ( const int vertexIndex : std::as_const( mInputVertices ) )
354 {
355 if ( meshEditor->isVertexFree( vertexIndex ) || meshEditor->isVertexOnBoundary( vertexIndex ) )
356 vertexIndextoTriangulate.append( vertexIndex );
357 else
358 removedVerticesFromTriangulation.append( vertexIndex );
359 }
360
361 bool triangulationReady = false;
362 bool giveUp = false;
364
365 while ( !triangulationReady )
366 {
367 QgsMeshTriangulation triangulation;
368
369 QVector<int> triangulationVertexToMeshVertex( vertexIndextoTriangulate.count() );
370 const QgsMesh *destinationMesh = meshEditor->topologicalMesh().mesh();
371
372 for ( int i = 0; i < vertexIndextoTriangulate.count(); ++i )
373 {
374 triangulationVertexToMeshVertex[i] = vertexIndextoTriangulate.at( i );
375 triangulation.addVertex( destinationMesh->vertices.at( vertexIndextoTriangulate.at( i ) ) );
376 }
377
378 QgsMesh resultingTriangulation = triangulation.triangulatedMesh();
379
380 //Transform the new mesh triangulation to destination mesh faces
381 QVector<QgsMeshFace> rawDestinationFaces = resultingTriangulation.faces;
382
383 for ( QgsMeshFace &destinationFace : rawDestinationFaces )
384 {
385 for ( int &vertexIndex : destinationFace )
386 vertexIndex = triangulationVertexToMeshVertex[vertexIndex];
387 }
388
389 //The new triangulation may contains faces that intersect existing faces, we need to remove them
390 QVector<QgsMeshFace> destinationFaces;
391 for ( const QgsMeshFace &face : rawDestinationFaces )
392 {
393 if ( meshEditor->isFaceGeometricallyCompatible( face ) )
394 destinationFaces.append( face );
395 }
396
397 bool facesReady = false;
398 QgsMeshEditingError previousError;
399 while ( !facesReady && !giveUp )
400 {
402 topologicFaces = QgsTopologicalMesh::createNewTopologicalFaces( destinationFaces, true, error );
403
404 if ( error == QgsMeshEditingError() )
405 error = meshEditor->topologicalMesh().facesCanBeAdded( topologicFaces );
406
407 switch ( error.errorType )
408 {
410 facesReady = true;
411 triangulationReady = true;
412 break;
417 if ( error.elementIndex != -1 )
418 destinationFaces.remove( error.elementIndex );
419 else
420 giveUp = true; //we don't know what happens, better to give up
421 break;
424 facesReady = true;
425 if ( error.elementIndex != -1 )
426 {
427 removedVerticesFromTriangulation.append( error.elementIndex );
428 vertexIndextoTriangulate.removeOne( error.elementIndex );
429 }
430 else
431 giveUp = true; //we don't know what happens, better to give up
432 break;
433 }
434 }
435 }
436
437 Q_ASSERT( meshEditor->topologicalMesh().checkConsistency() == QgsMeshEditingError() );
438
439 if ( !removedVerticesFromTriangulation.isEmpty() )
440 mMessage = QObject::tr( "%n vertices have not been included in the triangulation", nullptr, removedVerticesFromTriangulation.count() );
441
442 mIsFinished = true;
443
444 if ( triangulationReady && !giveUp )
445 return meshEditor->topologicalMesh().addFaces( topologicFaces );
446 else
448}
449
@ TooManyVerticesInFace
A face has more vertices than the maximum number supported per face.
@ InvalidFace
An error occurs due to an invalid face (for example, vertex indexes are unordered)
@ UniqueSharedVertex
A least two faces share only one vertices.
@ ManifoldFace
ManifoldFace.
@ InvalidVertex
An error occurs due to an invalid vertex (for example, vertex index is out of range the available ver...
@ FlatFace
A flat face is present.
GeometryType
The geometry types are used to group Qgis::WkbType in a coarse way.
Definition qgis.h:337
@ Polygon
Polygons.
@ PointZ
PointZ.
@ Forward
Forward transform (from source to destination)
The vertex_iterator class provides STL-style iterator for vertices.
This class represents a coordinate reference system (CRS).
Class for doing transforms between two map coordinate systems.
Custom exception class for Coordinate Reference System related exceptions.
Curve polygon geometry type.
Abstract base class for curved geometry type.
Definition qgscurve.h:35
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
QString what() const
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsGeometry geometry
Definition qgsfeature.h:69
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
A geometry is the spatial representation of a feature.
Qgis::GeometryOperationResult transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool transformZ=false)
Transforms this geometry as described by the coordinate transform ct.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Qgis::GeometryType type
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
@ SourceBreakLines
Break lines.
QgsMeshDataBlock is a block of integers/doubles that can be used to retrieve: active flags (e....
@ ScalarDouble
Scalar double values.
@ ActiveFlagInteger
Integer boolean flag whether face is active.
void setValues(const QVector< double > &vals)
Sets values.
void setValid(bool valid)
Sets block validity.
QgsMeshDatasetGroupMetadata is a collection of dataset group metadata such as whether the data is vec...
Abstract class that represents a dataset group.
void calculateStatistic() const
Calculates the statistics (minimum and maximum)
QgsMeshDatasetMetadata is a collection of mesh dataset metadata such as whether the data is valid or ...
QgsMeshDatasetValue represents single dataset value.
Abstract class that represents a dataset.
QString text() const override
Returns a short text string describing what this advanced edit does. Default implementation return a ...
Class that represents an error during mesh editing.
Qgis::MeshEditingErrorType errorType
Class that makes edit operation on a mesh.
bool isFaceGeometricallyCompatible(const QgsMeshFace &face) const
Returns true if the face does not intersect or contains any other elements (faces or vertices) The to...
bool isVertexOnBoundary(int vertexIndex) const
Returns whether the vertex with index vertexIndex is on a boundary.
bool isVertexFree(int vertexIndex) const
Returns whether the vertex with index vertexIndex is a free vertex.
QgsTopologicalMesh & topologicalMesh()
Returns a reference to the topological mesh.
Class that handles mesh creation with Delaunay constrained triangulation.
bool addBreakLines(QgsFeatureIterator &lineFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transformContext, QgsFeedback *feedback=nullptr, long featureCount=1)
Adds break lines from a vector layer, return true if successful.
int addVertex(const QgsPoint &vertex)
Adds a new vertex in the triangulation and returns the index of the new vertex.
void setCrs(const QgsCoordinateReferenceSystem &crs)
Sets the coordinate reference system used for the triangulation.
bool addVertices(QgsFeatureIterator &vertexFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback=nullptr, long featureCount=1)
Adds vertices to the triangulation from a feature iterator, return true if successful.
QgsMesh triangulatedMesh(QgsFeedback *feedback=nullptr) const
Returns the triangulated mesh.
int datasetCount() const override
Returns the count of datasets in the group.
QDomElement writeXml(QDomDocument &doc, const QgsReadWriteContext &context) const override
Write dataset group information in a DOM element.
void initialize() override
Initialize the dataset group.
QgsMeshDataset * dataset(int index) const override
Returns the dataset with index.
QgsMeshDatasetMetadata datasetMetadata(int datasetIndex) const override
Returns the metadata of the dataset with index datasetIndex.
QgsMeshZValueDatasetGroup(const QString &datasetGroupName, const QgsMesh &mesh)
Constructor.
QgsMeshDataBlock areFacesActive(int faceIndex, int count) const override
Returns whether faces are active.
QgsMeshDatasetMetadata metadata() const override
Returns the metadata of the dataset.
QgsMeshZValueDataset(const QgsMesh &mesh)
Constructor with the mesh.
bool isActive(int faceIndex) const override
Returns whether the face is active.
int valuesCount() const override
Returns the values count.
QgsMeshDatasetValue datasetValue(int valueIndex) const override
Returns the value with index valueIndex.
QgsMeshDataBlock datasetValues(bool isScalar, int valueIndex, int count) const override
Returns count values from valueIndex.
Multi curve geometry collection.
Multi surface geometry collection.
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
double y
Definition qgspoint.h:53
The class is used as a container of context for various read/write operations on other objects.
Class that contains topological differences between two states of a topological mesh,...
Class that contains independent faces an topological information about this faces.
QgsMeshEditingError checkConsistency() const
Checks the consistency of the topological mesh and return false if there is a consistency issue.
QgsMeshEditingError facesCanBeAdded(const TopologicalFaces &topologicalFaces) const
Returns whether the faces can be added to the mesh.
Changes addFaces(const TopologicalFaces &topologicFaces)
Adds faces topologicFaces to the topologic mesh.
QgsMesh * mesh() const
Returns a pointer to the wrapped mesh.
static TopologicalFaces createNewTopologicalFaces(const QVector< QgsMeshFace > &faces, bool uniqueSharedVertexAllowed, QgsMeshEditingError &error)
Creates new topological faces that are not yet included in the mesh.
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
QVector< QgsPoint > QgsPointSequence
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
QVector< int > QgsMeshFace
List of vertex indexes.
const QgsCoordinateReferenceSystem & crs
Mesh - vertices, edges and faces.
int vertexCount() const
Returns number of vertices.
QVector< QgsMeshVertex > vertices
QVector< QgsMeshFace > faces
int faceCount() const
Returns number of faces.
QgsMeshVertex vertex(int index) const
Returns a vertex at the index.