QGIS API Documentation  3.21.0-Master (5b68dc587e)
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 
18 #include "qgsmeshtriangulation.h"
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 
38 bool 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 
62 bool 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 
113 void 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 
147 void 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 
248 QgsMeshZValueDatasetGroup::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 
277 QDomElement 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 
303 QgsMeshDataBlock 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 
323 bool 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 
340 QgsTopologicalMesh::Changes QgsMeshEditingDelaunayTriangulation::apply( QgsMeshEditor *meshEditor )
341 {
342  //use only vertices that are on boundary or free, if boundary
343  QList<int> vertexIndextoTriangulate;
344 
345  QList<int> removedVerticesFromTriangulation;
346 
347  for ( const int vertexIndex : std::as_const( mInputVertices ) )
348  {
349  if ( meshEditor->isVertexFree( vertexIndex ) || meshEditor->isVertexOnBoundary( vertexIndex ) )
350  vertexIndextoTriangulate.append( vertexIndex );
351  else
352  removedVerticesFromTriangulation.append( vertexIndex );
353  }
354 
355  bool triangulationReady = false;
356  bool giveUp = false;
358 
359  while ( !triangulationReady )
360  {
361  QgsMeshTriangulation triangulation;
362 
363  QVector<int> triangulationVertexToMeshVertex( vertexIndextoTriangulate.count() );
364  const QgsMesh *destinationMesh = meshEditor->topologicalMesh().mesh();
365 
366  for ( int i = 0; i < vertexIndextoTriangulate.count(); ++i )
367  {
368  triangulationVertexToMeshVertex[i] = vertexIndextoTriangulate.at( i );
369  triangulation.addVertex( destinationMesh->vertices.at( vertexIndextoTriangulate.at( i ) ) );
370  }
371 
372  QgsMesh resultingTriangulation = triangulation.triangulatedMesh();
373 
374  //Transform the new mesh triangulation to destination mesh faces
375  QVector<QgsMeshFace> rawDestinationFaces = resultingTriangulation.faces;
376 
377  for ( QgsMeshFace &destinationFace : rawDestinationFaces )
378  {
379  for ( int &vertexIndex : destinationFace )
380  vertexIndex = triangulationVertexToMeshVertex[vertexIndex];
381  }
382 
383  //The new triangulation may contains faces that intersect existing faces, we need to remove them
384  QVector<QgsMeshFace> destinationFaces;
385  for ( const QgsMeshFace &face : rawDestinationFaces )
386  {
387  if ( meshEditor->isFaceGeometricallyCompatible( face ) )
388  destinationFaces.append( face );
389  }
390 
391  bool facesReady = false;
392  QgsMeshEditingError previousError;
393  while ( !facesReady && !giveUp )
394  {
395  QgsMeshEditingError error;
396  topologicFaces = meshEditor->topologicalMesh().createNewTopologicalFaces( destinationFaces, true, error );
397 
398  if ( error == QgsMeshEditingError() )
399  error = meshEditor->topologicalMesh().canFacesBeAdded( topologicFaces );
400 
401  switch ( error.errorType )
402  {
404  facesReady = true;
405  triangulationReady = true;
406  break;
411  if ( error.elementIndex != -1 )
412  destinationFaces.remove( error.elementIndex );
413  else
414  giveUp = true; //we don't know what happens, better to give up
415  break;
418  facesReady = true;
419  if ( error.elementIndex != -1 )
420  {
421  removedVerticesFromTriangulation.append( error.elementIndex );
422  vertexIndextoTriangulate.removeOne( error.elementIndex );
423  }
424  else
425  giveUp = true; //we don't know what happens, better to give up
426  break;
427  }
428  }
429  }
430 
431  Q_ASSERT( meshEditor->topologicalMesh().checkConsistency() == QgsMeshEditingError() );
432 
433  if ( !removedVerticesFromTriangulation.isEmpty() )
434  mMessage = QObject::tr( "%1 vertices have not been included in the triangulation" ).arg( removedVerticesFromTriangulation.count() );
435 
436  mIsFinished = true;
437 
438  if ( triangulationReady && !giveUp )
439  return meshEditor->topologicalMesh().addFaces( topologicFaces );
440  else
442 }
443 
@ 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:302
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:124
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:127
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()
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.
Class that represents an error during mesh editing.
Definition: qgsmesheditor.h:40
Qgis::MeshEditingErrorType errorType
Definition: qgsmesheditor.h:49
Class that makes edit operation on a mesh.
Definition: qgsmesheditor.h:65
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.
Changes addFaces(const TopologicalFaces &topologicFaces)
Adds faces topologicFaces to the topologic mesh.
QgsMesh * mesh() const
Returns a pointer to the wrapped mesh.
QgsMeshEditingError canFacesBeAdded(const TopologicalFaces &topologicalFaces) const
Returns whether the faces can be added to the 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:938
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.