QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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
80QgsGeometry 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
401void 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:164
QgsGeometry mergeLines() const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries, const QgsGeometryParameters &parameters=QgsGeometryParameters())
Compute the unary union on a list of geometries.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
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:100
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.
const QVector< QgsMeshFace > & triangles() const
Returns triangles.
const QVector< QgsMeshVertex > & vertices() const
Returns vertices in map coordinate system.
bool update(QgsMesh *nativeMesh, const QgsCoordinateTransform &transform)
Constructs triangular mesh from layer's native mesh and transform to destination CRS.
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:2527
QVector< int > QgsMeshFace
List of vertex indexes.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
QVector< QgsMeshFace > faces