QGIS API Documentation 3.29.0-Master (006c3c0232)
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
20#include "qgsvectorlayer.h"
22#include "qgscurve.h"
23#include "qgscurvepolygon.h"
24#include "qgsmultisurface.h"
25#include "qgsmulticurve.h"
26#include "qgsfeedback.h"
27#include "qgslogger.h"
28#include "qgsmesheditor.h"
29
31{
32 mTriangulation.reset( new QgsDualEdgeTriangulation() );
33}
34
35
37
38bool QgsMeshTriangulation::addVertices( QgsFeatureIterator &vertexFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback, long featureCount )
39{
40 if ( feedback )
41 feedback->setProgress( 0 );
42
43 QgsFeature feat;
44 long i = 0;
45 while ( vertexFeatureIterator.nextFeature( feat ) )
46 {
47 if ( feedback )
48 {
49 if ( feedback->isCanceled() )
50 break;
51
52 feedback->setProgress( 100 * i / featureCount );
53 i++;
54 }
55
56 addVerticesFromFeature( feat, valueAttribute, transform, feedback );
57 }
58
59 return true;
60}
61
62bool QgsMeshTriangulation::addBreakLines( QgsFeatureIterator &lineFeatureIterator, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback, long featureCount )
63{
64 if ( feedback )
65 feedback->setProgress( 0 );
66
67 QgsFeature feat;
68 long i = 0;
69 while ( lineFeatureIterator.nextFeature( feat ) )
70 {
71 if ( feedback )
72 {
73 if ( feedback->isCanceled() )
74 break;
75
76 feedback->setProgress( 100 * i / featureCount );
77 i++;
78 }
79
80 QgsWkbTypes::GeometryType geomType = feat.geometry().type();
81 switch ( geomType )
82 {
84 addVerticesFromFeature( feat, valueAttribute, transform, feedback );
85 break;
88 addBreakLinesFromFeature( feat, valueAttribute, transform, feedback );
89 break;
90 default:
91 break;
92 }
93 }
94
95 return true;
96}
97
99{
100 return mTriangulation->addPoint( vertex );
101}
102
104{
105 return mTriangulation->triangulationToMesh( feedback );
106}
107
109{
110 mCrs = crs;
111}
112
113void QgsMeshTriangulation::addVerticesFromFeature( const QgsFeature &feature, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback )
114{
115 QgsGeometry geom = feature.geometry();
116 try
117 {
118 geom.transform( transform, Qgis::TransformDirection::Forward, true );
119 }
120 catch ( QgsCsException &cse )
121 {
122 Q_UNUSED( cse )
123 QgsDebugMsgLevel( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ), 4 );
124 return;
125 }
126
128
129 double value = 0;
130 if ( valueAttribute >= 0 )
131 value = feature.attribute( valueAttribute ).toDouble();
132
133 while ( vit != geom.vertices_end() )
134 {
135 if ( feedback && feedback->isCanceled() )
136 break;
137 if ( valueAttribute < 0 )
138 mTriangulation->addPoint( *vit );
139 else
140 {
141 mTriangulation->addPoint( QgsPoint( QgsWkbTypes::PointZ, ( *vit ).x(), ( *vit ).y(), value ) );
142 }
143 ++vit;
144 }
145}
146
147void QgsMeshTriangulation::addBreakLinesFromFeature( const QgsFeature &feature, int valueAttribute, const QgsCoordinateTransform &transform, QgsFeedback *feedback )
148{
149 double valueOnVertex = 0;
150 if ( valueAttribute >= 0 )
151 valueOnVertex = feature.attribute( valueAttribute ).toDouble();
152
153 //from QgsTinInterpolator::insertData()
154 std::vector<const QgsCurve *> curves;
155 QgsGeometry geom = feature.geometry();
156 try
157 {
158 geom.transform( transform, Qgis::TransformDirection::Forward, true );
159 }
160 catch ( QgsCsException &cse )
161 {
162 Q_UNUSED( cse )
163 QgsDebugMsgLevel( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ), 4 );
164 return;
165 }
166
168 {
169 std::vector< const QgsCurvePolygon * > polygons;
170 if ( geom.isMultipart() )
171 {
172 const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( geom.constGet() );
173 for ( int i = 0; i < ms->numGeometries(); ++i )
174 {
175 polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( ms->geometryN( i ) ) );
176 }
177 }
178 else
179 {
180 polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( geom.constGet() ) );
181 }
182
183 for ( const QgsCurvePolygon *polygon : polygons )
184 {
185 if ( feedback && feedback->isCanceled() )
186 break;
187 if ( !polygon )
188 continue;
189
190 if ( polygon->exteriorRing() )
191 curves.emplace_back( polygon->exteriorRing() );
192
193 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
194 {
195 if ( feedback && feedback->isCanceled() )
196 break;
197 curves.emplace_back( polygon->interiorRing( i ) );
198 }
199 }
200 }
201 else
202 {
203 if ( geom.isMultipart() )
204 {
205 const QgsMultiCurve *mc = qgsgeometry_cast< const QgsMultiCurve * >( geom.constGet() );
206 for ( int i = 0; i < mc->numGeometries(); ++i )
207 {
208 if ( feedback && feedback->isCanceled() )
209 break;
210 curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( mc->geometryN( i ) ) );
211 }
212 }
213 else
214 {
215 curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( geom.constGet() ) );
216 }
217 }
218
219 for ( const QgsCurve *curve : curves )
220 {
221 if ( !curve )
222 continue;
223
224 if ( feedback && feedback->isCanceled() )
225 break;
226
227 QgsPointSequence linePoints;
228 curve->points( linePoints );
229 bool hasZ = curve->is3D();
230 if ( valueAttribute >= 0 )
231 for ( int i = 0; i < linePoints.count(); ++i )
232 {
233 if ( feedback && feedback->isCanceled() )
234 break;
235 if ( hasZ )
236 linePoints[i].setZ( valueOnVertex );
237 else
238 {
239 const QgsPoint &point = linePoints.at( i );
240 linePoints[i] = QgsPoint( point.x(), point.y(), valueOnVertex );
241 }
242 }
243
244 mTriangulation->addLine( linePoints, QgsInterpolator::SourceBreakLines );
245 }
246}
247
248QgsMeshZValueDatasetGroup::QgsMeshZValueDatasetGroup( const QString &datasetGroupName, const QgsMesh &mesh ):
249 QgsMeshDatasetGroup( datasetGroupName, QgsMeshDatasetGroupMetadata::DataOnVertices )
250{
251 mDataset = std::make_unique< QgsMeshZValueDataset >( mesh );
252}
253
255{
257}
258
260{
261 if ( datasetIndex != 0 )
262 return QgsMeshDatasetMetadata();
263
264 return mDataset->metadata();
265}
266
268
270{
271 if ( index != 0 )
272 return nullptr;
273
274 return mDataset.get();
275}
276
277QDomElement QgsMeshZValueDatasetGroup::writeXml( QDomDocument &doc, const QgsReadWriteContext &context ) const
278{
279 Q_UNUSED( doc );
280 Q_UNUSED( context );
281 return QDomElement();
282}
283
285{
286 for ( const QgsMeshVertex &vertex : mesh.vertices )
287 {
288 if ( vertex.z() < mZMinimum )
289 mZMinimum = vertex.z();
290 if ( vertex.z() > mZMaximum )
291 mZMaximum = vertex.z();
292 }
293}
294
296{
297 if ( valueIndex < 0 || valueIndex >= mMesh.vertexCount() )
298 return QgsMeshDatasetValue();
299
300 return QgsMeshDatasetValue( mMesh.vertex( valueIndex ).z() );
301}
302
303QgsMeshDataBlock QgsMeshZValueDataset::datasetValues( bool isScalar, int valueIndex, int count ) const
304{
305 Q_UNUSED( isScalar )
307 QVector<double> zValues( count );
308 for ( int i = valueIndex; i < valueIndex + count; ++i )
309 zValues[i - valueIndex] = mMesh.vertex( i ).z();
310 ret.setValues( zValues );
311 return ret;
312}
313
315{
316 Q_UNUSED( faceIndex );
317 Q_UNUSED( count );
319 ret.setValid( true );
320 return ret;
321}
322
323bool QgsMeshZValueDataset::isActive( int faceIndex ) const
324{
325 return ( faceIndex > 0 && faceIndex < mMesh.faceCount() );
326}
327
329{
330 return QgsMeshDatasetMetadata( 0, true, mZMinimum, mZMaximum, 0 );
331}
332
334{
335 return mMesh.vertexCount();
336}
337
339
341{
342 return QObject::tr( "Delaunay triangulation" );
343}
344
345QgsTopologicalMesh::Changes QgsMeshEditingDelaunayTriangulation::apply( QgsMeshEditor *meshEditor )
346{
347 //use only vertices that are on boundary or free, if boundary
348 QList<int> vertexIndextoTriangulate;
349
350 QList<int> removedVerticesFromTriangulation;
351
352 for ( const int vertexIndex : std::as_const( mInputVertices ) )
353 {
354 if ( meshEditor->isVertexFree( vertexIndex ) || meshEditor->isVertexOnBoundary( vertexIndex ) )
355 vertexIndextoTriangulate.append( vertexIndex );
356 else
357 removedVerticesFromTriangulation.append( vertexIndex );
358 }
359
360 bool triangulationReady = false;
361 bool giveUp = false;
363
364 while ( !triangulationReady )
365 {
366 QgsMeshTriangulation triangulation;
367
368 QVector<int> triangulationVertexToMeshVertex( vertexIndextoTriangulate.count() );
369 const QgsMesh *destinationMesh = meshEditor->topologicalMesh().mesh();
370
371 for ( int i = 0; i < vertexIndextoTriangulate.count(); ++i )
372 {
373 triangulationVertexToMeshVertex[i] = vertexIndextoTriangulate.at( i );
374 triangulation.addVertex( destinationMesh->vertices.at( vertexIndextoTriangulate.at( i ) ) );
375 }
376
377 QgsMesh resultingTriangulation = triangulation.triangulatedMesh();
378
379 //Transform the new mesh triangulation to destination mesh faces
380 QVector<QgsMeshFace> rawDestinationFaces = resultingTriangulation.faces;
381
382 for ( QgsMeshFace &destinationFace : rawDestinationFaces )
383 {
384 for ( int &vertexIndex : destinationFace )
385 vertexIndex = triangulationVertexToMeshVertex[vertexIndex];
386 }
387
388 //The new triangulation may contains faces that intersect existing faces, we need to remove them
389 QVector<QgsMeshFace> destinationFaces;
390 for ( const QgsMeshFace &face : rawDestinationFaces )
391 {
392 if ( meshEditor->isFaceGeometricallyCompatible( face ) )
393 destinationFaces.append( face );
394 }
395
396 bool facesReady = false;
397 QgsMeshEditingError previousError;
398 while ( !facesReady && !giveUp )
399 {
401 topologicFaces = QgsTopologicalMesh::createNewTopologicalFaces( destinationFaces, true, error );
402
403 if ( error == QgsMeshEditingError() )
404 error = meshEditor->topologicalMesh().facesCanBeAdded( topologicFaces );
405
406 switch ( error.errorType )
407 {
409 facesReady = true;
410 triangulationReady = true;
411 break;
416 if ( error.elementIndex != -1 )
417 destinationFaces.remove( error.elementIndex );
418 else
419 giveUp = true; //we don't know what happens, better to give up
420 break;
423 facesReady = true;
424 if ( error.elementIndex != -1 )
425 {
426 removedVerticesFromTriangulation.append( error.elementIndex );
427 vertexIndextoTriangulate.removeOne( error.elementIndex );
428 }
429 else
430 giveUp = true; //we don't know what happens, better to give up
431 break;
432 }
433 }
434 }
435
436 Q_ASSERT( meshEditor->topologicalMesh().checkConsistency() == QgsMeshEditingError() );
437
438 if ( !removedVerticesFromTriangulation.isEmpty() )
439 mMessage = QObject::tr( "%n vertices have not been included in the triangulation", nullptr, removedVerticesFromTriangulation.count() );
440
441 mIsFinished = true;
442
443 if ( triangulationReady && !giveUp )
444 return meshEditor->topologicalMesh().addFaces( topologicFaces );
445 else
447}
448
@ 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.
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.
Definition: qgsexception.h:66
Curve polygon geometry type.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
QString what() const
Definition: qgsexception.h:48
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsGeometry geometry
Definition: qgsfeature.h:67
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:338
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
int numGeometries() const SIP_HOLDGIL
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.
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Qgis::GeometryOperationResult transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool transformZ=false) SIP_THROW(QgsCsException)
Transforms this geometry as described by the coordinate transform ct.
bool isMultipart() const SIP_HOLDGIL
Returns true if WKB of the geometry is of WKBMulti* type.
QgsWkbTypes::GeometryType type
Definition: qgsgeometry.h:167
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
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.
QgsMeshEditingDelaunayTriangulation()
Constructor.
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.
Definition: qgsmesheditor.h:43
Qgis::MeshEditingErrorType errorType
Definition: qgsmesheditor.h:52
Class that makes edit operation on a mesh.
Definition: qgsmesheditor.h:68
bool isVertexOnBoundary(int vertexIndex) const
Returns whether the vertex with index vertexIndex is on a boundary.
bool isFaceGeometricallyCompatible(const QgsMeshFace &face)
Returns true if the face does not intersect or contains any other elements (faces or vertices) The to...
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.
~QgsMeshTriangulation()
Destructor.
QgsMeshTriangulation()
Constructor.
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.
Definition: qgsmulticurve.h:30
Multi surface geometry collection.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Q_GADGET double x
Definition: qgspoint.h:52
double z
Definition: qgspoint.h:54
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 GeometryType geometryType(Type type) SIP_HOLDGIL
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:968
GeometryType
The geometry types are used to group QgsWkbTypes::Type in a coarse way.
Definition: qgswkbtypes.h:141
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.