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