QGIS API Documentation  3.23.0-Master (dd0cd13a00)
qgsmeshcontours.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsmeshcontours.cpp
3  -------------------
4  begin : September 25th, 2019
5  copyright : (C) 2018 by Peter Petrik
6  email : zilolv 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 #include "qgsmeshcontours.h"
18 #include "qgspoint.h"
19 #include "qgsmultilinestring.h"
20 #include "qgslinestring.h"
21 #include "qgspolygon.h"
22 #include "qgsmapsettings.h"
23 #include "qgsmeshlayer.h"
24 #include "qgstriangularmesh.h"
25 #include "qgsgeometryutils.h"
26 #include "qgsfeature.h"
27 #include "qgsgeometry.h"
28 #include "qgsmeshlayerutils.h"
29 #include "qgsmeshdataprovider.h"
30 #include "qgsfeedback.h"
31 #include "qgsproject.h"
32 
33 #include <limits>
34 
35 #include <QSet>
36 #include <QPair>
37 
39  : mMeshLayer( layer )
40 {
41  if ( !mMeshLayer || !mMeshLayer->dataProvider() || !mMeshLayer->dataProvider()->isValid() )
42  return;
43 
44  // Support for meshes with edges is not implemented
45  if ( mMeshLayer->dataProvider()->contains( QgsMesh::ElementType::Edge ) )
46  return;
47 
48  mMeshLayer->dataProvider()->populateMesh( &mNativeMesh );
49  mTriangularMesh.update( &mNativeMesh );
50 }
51 
53  const QgsMesh &nativeMesh,
54  const QVector<double> &datasetValues,
55  const QgsMeshDataBlock scalarActiveFaceFlagValues ):
56  mTriangularMesh( triangularMesh )
57  , mNativeMesh( nativeMesh )
58  , mDatasetValues( datasetValues )
59  , mScalarActiveFaceFlagValues( scalarActiveFaceFlagValues )
60 {}
61 
63 
65  const QgsMeshDatasetIndex &index,
66  double min_value,
67  double max_value,
69  QgsFeedback *feedback
70 )
71 {
72  if ( !mMeshLayer )
73  return QgsGeometry();
74 
75  populateCache( index, method );
76 
77  return exportPolygons( min_value, max_value, feedback );
78 }
79 
80 QgsGeometry QgsMeshContours::exportPolygons( double min_value, double max_value, QgsFeedback *feedback )
81 {
82  if ( min_value > max_value )
83  {
84  const double tmp = max_value;
85  max_value = min_value;
86  min_value = tmp;
87  }
88 
89  QVector<QgsGeometry> multiPolygon;
90 
91  // STEP 1: Get Data
92  const QVector<QgsMeshVertex> vertices = mTriangularMesh.vertices();
93  const QVector<int> &trianglesToNativeFaces = mTriangularMesh.trianglesToNativeFaces();
94 
95  // STEP 2: For each triangle get the contour line
96  for ( int i = 0; i < mTriangularMesh.triangles().size(); ++i )
97  {
98  if ( feedback && feedback->isCanceled() )
99  break;
100 
101  const int nativeIndex = trianglesToNativeFaces.at( i );
102  if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
103  continue;
104 
105  const QgsMeshFace &triangle = mTriangularMesh.triangles().at( i );
106  const int indices[3] =
107  {
108  triangle.at( 0 ),
109  triangle.at( 1 ),
110  triangle.at( 2 )
111  };
112 
113  const QVector<QgsMeshVertex> coords =
114  {
115  vertices.at( indices[0] ),
116  vertices.at( indices[1] ),
117  vertices.at( indices[2] )
118  };
119 
120  const double values[3] =
121  {
122  mDatasetValues.at( indices[0] ),
123  mDatasetValues.at( indices[1] ),
124  mDatasetValues.at( indices[2] )
125  };
126 
127  // any value is NaN
128  if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
129  continue;
130 
131  // all values on vertices are outside the range
132  if ( ( ( min_value > values[0] ) && ( min_value > values[1] ) && ( min_value > values[2] ) ) ||
133  ( ( max_value < values[0] ) && ( max_value < values[1] ) && ( max_value < values[2] ) ) )
134  continue;
135 
136  const bool valueInRange[3] =
137  {
138  ( min_value <= values[0] ) &&( max_value >= values[0] ),
139  ( min_value <= values[1] ) &&( max_value >= values[1] ),
140  ( min_value <= values[2] ) &&( max_value >= values[2] )
141  };
142 
143  // all values are inside the range == take whole triangle
144  if ( valueInRange[0] && valueInRange[1] && valueInRange[2] )
145  {
146  QVector<QgsMeshVertex> ring = coords;
147  ring.push_back( coords[0] );
148  std::unique_ptr< QgsLineString > ext = std::make_unique< QgsLineString> ( coords );
149  std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
150  poly->setExteriorRing( ext.release() );
151  multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
152  continue;
153  }
154 
155  // go through all edges
156  QVector<QgsMeshVertex> ring;
157  for ( int i = 0; i < 3; ++i )
158  {
159  const int j = ( i + 1 ) % 3;
160 
161  if ( valueInRange[i] )
162  {
163  if ( valueInRange[j] )
164  {
165  // whole edge is part of resulting contour polygon edge
166  if ( !ring.contains( coords[i] ) )
167  ring.push_back( coords[i] );
168  if ( !ring.contains( coords[j] ) )
169  ring.push_back( coords[j] );
170  }
171  else
172  {
173  // i is part or the resulting edge
174  if ( !ring.contains( coords[i] ) )
175  ring.push_back( coords[i] );
176  // we need to find the other point
177  double value = max_value;
178  if ( values[i] > values[j] )
179  {
180  value = min_value;
181  }
182  const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
183  const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
184  if ( !ring.contains( xy ) )
185  ring.push_back( xy );
186  }
187  }
188  else
189  {
190  if ( valueInRange[j] )
191  {
192  // we need to find the other point
193  double value = max_value;
194  if ( values[i] < values[j] )
195  {
196  value = min_value;
197  }
198 
199  const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
200  const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
201  if ( !ring.contains( xy ) )
202  ring.push_back( xy );
203 
204  // j is part
205  if ( !ring.contains( coords[j] ) )
206  ring.push_back( coords[j] );
207 
208  }
209  else
210  {
211  // last option we need to consider is that both min and max are between
212  // value i and j, and in that case we need to calculate both point
213  double value1 = max_value;
214  double value2 = max_value;
215  if ( values[i] < values[j] )
216  {
217  if ( ( min_value < values[i] ) || ( max_value > values[j] ) )
218  {
219  continue;
220  }
221  value1 = min_value;
222  }
223  else
224  {
225  if ( ( min_value < values[j] ) || ( max_value > values[i] ) )
226  {
227  continue;
228  }
229  value2 = min_value;
230  }
231 
232  const double fraction1 = ( value1 - values[i] ) / ( values[j] - values[i] );
233  const QgsPoint xy1 = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction1 );
234  if ( !ring.contains( xy1 ) )
235  ring.push_back( xy1 );
236 
237  const double fraction2 = ( value2 - values[i] ) / ( values[j] - values[i] );
238  const QgsPoint xy2 = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction2 );
239  if ( !ring.contains( xy2 ) )
240  ring.push_back( xy2 );
241  }
242  }
243  }
244 
245  // add if the polygon is not degraded
246  if ( ring.size() > 2 )
247  {
248  std::unique_ptr< QgsLineString > ext = std::make_unique< QgsLineString> ( ring );
249  std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
250  poly->setExteriorRing( ext.release() );
251  multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
252  }
253  }
254 
255  // STEP 3: dissolve the individual polygons from triangles if possible
256  if ( multiPolygon.isEmpty() )
257  {
258  return QgsGeometry();
259  }
260  else
261  {
262  const QgsGeometry res = QgsGeometry::unaryUnion( multiPolygon );
263  return res;
264  }
265 }
266 
268  double value,
270  QgsFeedback *feedback )
271 {
272  // Check if the layer/mesh is valid
273  if ( !mMeshLayer )
274  return QgsGeometry();
275 
276  populateCache( index, method );
277 
278  return exportLines( value, feedback );
279 }
280 
282 {
283  std::unique_ptr<QgsMultiLineString> multiLineString( new QgsMultiLineString() );
284  QSet<QPair<int, int>> exactEdges;
285 
286  // STEP 1: Get Data
287  const QVector<QgsMeshVertex> vertices = mTriangularMesh.vertices();
288  const QVector<int> &trianglesToNativeFaces = mTriangularMesh.trianglesToNativeFaces();
289 
290  // STEP 2: For each triangle get the contour line
291  for ( int i = 0; i < mTriangularMesh.triangles().size(); ++i )
292  {
293  if ( feedback && feedback->isCanceled() )
294  break;
295 
296  const int nativeIndex = trianglesToNativeFaces.at( i );
297  if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
298  continue;
299 
300  const QgsMeshFace &triangle = mTriangularMesh.triangles().at( i );
301 
302  const int indices[3] =
303  {
304  triangle.at( 0 ),
305  triangle.at( 1 ),
306  triangle.at( 2 )
307  };
308 
309  const QVector<QgsMeshVertex> coords =
310  {
311  vertices.at( indices[0] ),
312  vertices.at( indices[1] ),
313  vertices.at( indices[2] )
314  };
315 
316  const double values[3] =
317  {
318  mDatasetValues.at( indices[0] ),
319  mDatasetValues.at( indices[1] ),
320  mDatasetValues.at( indices[2] )
321  };
322 
323  // any value is NaN
324  if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
325  continue;
326 
327  // value is outside the range
328  if ( ( ( value > values[0] ) && ( value > values[1] ) && ( value > values[2] ) ) ||
329  ( ( value < values[0] ) && ( value < values[1] ) && ( value < values[2] ) ) )
330  continue;
331 
332  // all values are the same
333  if ( qgsDoubleNear( values[0], values[1] ) && qgsDoubleNear( values[1], values[2] ) )
334  continue;
335 
336  // go through all edges
337  QgsPoint tmp;
338 
339  for ( int i = 0; i < 3; ++i )
340  {
341  const int j = ( i + 1 ) % 3;
342  // value is outside the range
343  if ( ( ( value > values[i] ) && ( value > values[j] ) ) ||
344  ( ( value < values[i] ) && ( value < values[j] ) ) )
345  continue;
346 
347  // the whole edge is result and we are done
348  if ( qgsDoubleNear( values[i], values[j] ) && qgsDoubleNear( values[i], values[j] ) )
349  {
350  if ( exactEdges.contains( { indices[i], indices[j] } ) || exactEdges.contains( { indices[j], indices[i] } ) )
351  {
352  break;
353  }
354  else
355  {
356  exactEdges.insert( { indices[i], indices[j] } );
357  std::unique_ptr<QgsLineString> line( new QgsLineString( coords[i], coords[j] ) );
358  multiLineString->addGeometry( line.release() );
359  break;
360  }
361  }
362 
363  // only one point matches, we are not interested in this
364  if ( qgsDoubleNear( values[i], value ) || qgsDoubleNear( values[j], value ) )
365  {
366  continue;
367  }
368 
369  // ok part of the result contour line is one point on this edge
370  const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
371  const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
372 
373  if ( std::isnan( tmp.x() ) )
374  {
375  // ok we have found start point of the contour line
376  tmp = xy;
377  }
378  else
379  {
380  // we have found the end point of the contour line, we are done
381  std::unique_ptr<QgsLineString> line( new QgsLineString( tmp, xy ) );
382  multiLineString->addGeometry( line.release() );
383  break;
384  }
385  }
386  }
387 
388  // STEP 3: merge the contour segments to linestrings
389  if ( multiLineString->isEmpty() )
390  {
391  return QgsGeometry();
392  }
393  else
394  {
395  const QgsGeometry in( multiLineString.release() );
396  const QgsGeometry res = in.mergeLines();
397  return res;
398  }
399 }
400 
401 void QgsMeshContours::populateCache( const QgsMeshDatasetIndex &index, QgsMeshRendererScalarSettings::DataResamplingMethod method )
402 {
403  if ( !mMeshLayer )
404  return;
405 
406  if ( mCachedIndex != index )
407  {
408  const bool scalarDataOnVertices = mMeshLayer->dataProvider()->datasetGroupMetadata( index ).dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
409  const int count = scalarDataOnVertices ? mNativeMesh.vertices.count() : mNativeMesh.faces.count();
410 
411  // populate scalar values
412  const QgsMeshDataBlock vals = QgsMeshLayerUtils::datasetValues(
413  mMeshLayer,
414  index,
415  0,
416  count );
417  if ( vals.isValid() )
418  {
419  // vals could be scalar or vectors, for contour rendering we want always magnitude
420  mDatasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
421  }
422  else
423  {
424  mDatasetValues = QVector<double>( count, std::numeric_limits<double>::quiet_NaN() );
425  }
426 
427  // populate face active flag, always defined on faces
428  mScalarActiveFaceFlagValues = mMeshLayer->dataProvider()->areFacesActive(
429  index,
430  0,
431  mNativeMesh.faces.count() );
432 
433  // for data on faces, there could be request to interpolate the data to vertices
434  if ( ( !scalarDataOnVertices ) )
435  {
436  mDatasetValues = QgsMeshLayerUtils::interpolateFromFacesData(
437  mDatasetValues,
438  &mNativeMesh,
439  &mTriangularMesh,
440  &mScalarActiveFaceFlagValues,
441  method
442  );
443  }
444  }
445 }
virtual bool isValid() const =0
Returns true if this is a valid layer.
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
static QgsPointXY interpolatePointOnLine(double x1, double y1, double x2, double y2, double fraction) SIP_HOLDGIL
Interpolates the position of a point a fraction of the way along the line from (x1,...
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:125
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries)
Compute the unary union on a list of geometries.
QgsGeometry mergeLines() const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:44
QgsGeometry exportPolygons(const QgsMeshDatasetIndex &index, double min_value, double max_value, QgsMeshRendererScalarSettings::DataResamplingMethod method, QgsFeedback *feedback=nullptr)
Exports multi polygons representing the areas with values in range for particular dataset.
QgsGeometry exportLines(const QgsMeshDatasetIndex &index, double value, QgsMeshRendererScalarSettings::DataResamplingMethod method, QgsFeedback *feedback=nullptr)
Exports multi line string containing the contour line for particular dataset and value.
QgsMeshContours(QgsMeshLayer *layer)
Constructs the mesh contours exporter.
QgsMeshDataBlock is a block of integers/doubles that can be used to retrieve: active flags (e....
bool isValid() const
Whether the block is valid.
bool active(int index) const
Returns a value for active flag by the index For scalar and vector 2d the behavior is undefined.
bool contains(const QgsMesh::ElementType &type) const
Returns whether the mesh contains at mesh elements of given type.
virtual void populateMesh(QgsMesh *mesh) const =0
Populates the mesh vertices, edges and faces.
DataType dataType() const
Returns whether dataset group data is defined on vertices or faces or volumes.
@ DataOnVertices
Data is defined on vertices.
QgsMeshDatasetIndex is index that identifies the dataset group (e.g.
virtual QgsMeshDatasetGroupMetadata datasetGroupMetadata(int groupIndex) const =0
Returns dataset group metadata.
virtual QgsMeshDataBlock areFacesActive(QgsMeshDatasetIndex index, int faceIndex, int count) const =0
Returns whether the faces are active for particular dataset.
Represents a mesh layer supporting display of data on structured or unstructured meshes.
Definition: qgsmeshlayer.h:97
QgsMeshDataProvider * dataProvider() override
Returns the layer's data provider, it may be nullptr.
DataResamplingMethod
Resampling of value from dataset.
Multi line string 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
Triangular/Derived Mesh is mesh with vertices in map coordinates.
bool update(QgsMesh *nativeMesh, const QgsCoordinateTransform &transform=QgsCoordinateTransform())
Constructs triangular mesh from layer's native mesh and transform to destination CRS.
const QVector< QgsMeshFace > & triangles() const
Returns triangles.
const QVector< QgsMeshVertex > & vertices() const
Returns vertices in map coordinate system.
const QVector< int > & trianglesToNativeFaces() const
Returns mapping between triangles and original faces.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1504
QVector< int > QgsMeshFace
List of vertex indexes.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
QVector< QgsMeshFace > faces