Branch data Line data Source code
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 : :
38 : 0 : QgsMeshContours::QgsMeshContours( QgsMeshLayer *layer )
39 : 0 : : mMeshLayer( layer )
40 : : {
41 : 0 : if ( !mMeshLayer || !mMeshLayer->dataProvider() || !mMeshLayer->dataProvider()->isValid() )
42 : 0 : return;
43 : :
44 : : // Support for meshes with edges is not implemented
45 : 0 : if ( mMeshLayer->dataProvider()->contains( QgsMesh::ElementType::Edge ) )
46 : 0 : return;
47 : :
48 : 0 : mMeshLayer->dataProvider()->populateMesh( &mNativeMesh );
49 : 0 : mTriangularMesh.update( &mNativeMesh );
50 : 0 : }
51 : :
52 : 0 : QgsMeshContours::QgsMeshContours( const QgsTriangularMesh &triangularMesh,
53 : : const QgsMesh &nativeMesh,
54 : : const QVector<double> &datasetValues,
55 : : const QgsMeshDataBlock scalarActiveFaceFlagValues ):
56 : 0 : mTriangularMesh( triangularMesh )
57 : 0 : , mNativeMesh( nativeMesh )
58 : 0 : , mDatasetValues( datasetValues )
59 : 0 : , mScalarActiveFaceFlagValues( scalarActiveFaceFlagValues )
60 : 0 : {}
61 : :
62 : 0 : QgsMeshContours::~QgsMeshContours() = default;
63 : :
64 : 0 : QgsGeometry QgsMeshContours::exportPolygons(
65 : : const QgsMeshDatasetIndex &index,
66 : : double min_value,
67 : : double max_value,
68 : : QgsMeshRendererScalarSettings::DataResamplingMethod method,
69 : : QgsFeedback *feedback
70 : : )
71 : : {
72 : 0 : if ( !mMeshLayer )
73 : 0 : return QgsGeometry();
74 : :
75 : 0 : populateCache( index, method );
76 : :
77 : 0 : return exportPolygons( min_value, max_value, feedback );
78 : 0 : }
79 : :
80 : 0 : QgsGeometry QgsMeshContours::exportPolygons( double min_value, double max_value, QgsFeedback *feedback )
81 : : {
82 : 0 : if ( min_value > max_value )
83 : : {
84 : 0 : double tmp = max_value;
85 : 0 : max_value = min_value;
86 : 0 : min_value = tmp;
87 : 0 : }
88 : :
89 : 0 : QVector<QgsGeometry> multiPolygon;
90 : :
91 : : // STEP 1: Get Data
92 : 0 : const QVector<QgsMeshVertex> vertices = mTriangularMesh.vertices();
93 : 0 : const QVector<int> &trianglesToNativeFaces = mTriangularMesh.trianglesToNativeFaces();
94 : :
95 : : // STEP 2: For each triangle get the contour line
96 : 0 : for ( int i = 0; i < mTriangularMesh.triangles().size(); ++i )
97 : : {
98 : 0 : if ( feedback && feedback->isCanceled() )
99 : 0 : break;
100 : :
101 : 0 : int nativeIndex = trianglesToNativeFaces.at( i );
102 : 0 : if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
103 : 0 : continue;
104 : :
105 : 0 : const QgsMeshFace &triangle = mTriangularMesh.triangles().at( i );
106 : : const int indices[3] =
107 : 0 : {
108 : 0 : triangle.at( 0 ),
109 : 0 : triangle.at( 1 ),
110 : 0 : triangle.at( 2 )
111 : : };
112 : :
113 : : const QVector<QgsMeshVertex> coords =
114 : 0 : {
115 : 0 : vertices.at( indices[0] ),
116 : 0 : vertices.at( indices[1] ),
117 : 0 : vertices.at( indices[2] )
118 : : };
119 : :
120 : : const double values[3] =
121 : 0 : {
122 : 0 : mDatasetValues.at( indices[0] ),
123 : 0 : mDatasetValues.at( indices[1] ),
124 : 0 : mDatasetValues.at( indices[2] )
125 : : };
126 : :
127 : : // any value is NaN
128 : 0 : if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
129 : 0 : continue;
130 : :
131 : : // all values on vertices are outside the range
132 : 0 : if ( ( ( min_value > values[0] ) && ( min_value > values[1] ) && ( min_value > values[2] ) ) ||
133 : 0 : ( ( max_value < values[0] ) && ( max_value < values[1] ) && ( max_value < values[2] ) ) )
134 : 0 : continue;
135 : :
136 : : const bool valueInRange[3] =
137 : 0 : {
138 : 0 : ( min_value <= values[0] ) &&( max_value >= values[0] ),
139 : 0 : ( min_value <= values[1] ) &&( max_value >= values[1] ),
140 : 0 : ( min_value <= values[2] ) &&( max_value >= values[2] )
141 : : };
142 : :
143 : : // all values are inside the range == take whole triangle
144 : 0 : if ( valueInRange[0] && valueInRange[1] && valueInRange[2] )
145 : 0 : {
146 : 0 : QVector<QgsMeshVertex> ring = coords;
147 : 0 : ring.push_back( coords[0] );
148 : 0 : std::unique_ptr< QgsLineString > ext = std::make_unique< QgsLineString> ( coords );
149 : 0 : std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
150 : 0 : poly->setExteriorRing( ext.release() );
151 : 0 : multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
152 : : continue;
153 : 0 : }
154 : :
155 : : // go through all edges
156 : 0 : QVector<QgsMeshVertex> ring;
157 : 0 : for ( int i = 0; i < 3; ++i )
158 : : {
159 : 0 : const int j = ( i + 1 ) % 3;
160 : :
161 : 0 : if ( valueInRange[i] )
162 : : {
163 : 0 : if ( valueInRange[j] )
164 : : {
165 : : // whole edge is part of resulting contour polygon edge
166 : 0 : if ( !ring.contains( coords[i] ) )
167 : 0 : ring.push_back( coords[i] );
168 : 0 : if ( !ring.contains( coords[j] ) )
169 : 0 : ring.push_back( coords[j] );
170 : 0 : }
171 : : else
172 : : {
173 : : // i is part or the resulting edge
174 : 0 : if ( !ring.contains( coords[i] ) )
175 : 0 : ring.push_back( coords[i] );
176 : : // we need to find the other point
177 : 0 : double value = max_value;
178 : 0 : if ( values[i] > values[j] )
179 : : {
180 : 0 : value = min_value;
181 : 0 : }
182 : 0 : const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
183 : 0 : const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
184 : 0 : if ( !ring.contains( xy ) )
185 : 0 : ring.push_back( xy );
186 : 0 : }
187 : 0 : }
188 : : else
189 : : {
190 : 0 : if ( valueInRange[j] )
191 : : {
192 : : // we need to find the other point
193 : 0 : double value = max_value;
194 : 0 : if ( values[i] < values[j] )
195 : : {
196 : 0 : value = min_value;
197 : 0 : }
198 : :
199 : 0 : const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
200 : 0 : const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
201 : 0 : if ( !ring.contains( xy ) )
202 : 0 : ring.push_back( xy );
203 : :
204 : : // j is part
205 : 0 : if ( !ring.contains( coords[j] ) )
206 : 0 : ring.push_back( coords[j] );
207 : :
208 : 0 : }
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 : 0 : double value1 = max_value;
214 : 0 : double value2 = max_value;
215 : 0 : if ( values[i] < values[j] )
216 : : {
217 : 0 : if ( ( min_value < values[i] ) || ( max_value > values[j] ) )
218 : : {
219 : 0 : continue;
220 : : }
221 : 0 : value1 = min_value;
222 : 0 : }
223 : : else
224 : : {
225 : 0 : if ( ( min_value < values[j] ) || ( max_value > values[i] ) )
226 : : {
227 : 0 : continue;
228 : : }
229 : 0 : value2 = min_value;
230 : : }
231 : :
232 : 0 : const double fraction1 = ( value1 - values[i] ) / ( values[j] - values[i] );
233 : 0 : const QgsPoint xy1 = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction1 );
234 : 0 : if ( !ring.contains( xy1 ) )
235 : 0 : ring.push_back( xy1 );
236 : :
237 : 0 : const double fraction2 = ( value2 - values[i] ) / ( values[j] - values[i] );
238 : 0 : const QgsPoint xy2 = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction2 );
239 : 0 : if ( !ring.contains( xy2 ) )
240 : 0 : ring.push_back( xy2 );
241 : 0 : }
242 : : }
243 : 0 : }
244 : :
245 : : // add if the polygon is not degraded
246 : 0 : if ( ring.size() > 2 )
247 : : {
248 : 0 : std::unique_ptr< QgsLineString > ext = std::make_unique< QgsLineString> ( ring );
249 : 0 : std::unique_ptr< QgsPolygon > poly = std::make_unique< QgsPolygon >();
250 : 0 : poly->setExteriorRing( ext.release() );
251 : 0 : multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
252 : 0 : }
253 : 0 : }
254 : :
255 : : // STEP 3: dissolve the individual polygons from triangles if possible
256 : 0 : if ( multiPolygon.isEmpty() )
257 : : {
258 : 0 : return QgsGeometry();
259 : : }
260 : : else
261 : : {
262 : 0 : QgsGeometry res = QgsGeometry::unaryUnion( multiPolygon );
263 : 0 : return res;
264 : 0 : }
265 : 0 : }
266 : :
267 : 0 : QgsGeometry QgsMeshContours::exportLines( const QgsMeshDatasetIndex &index,
268 : : double value,
269 : : QgsMeshRendererScalarSettings::DataResamplingMethod method,
270 : : QgsFeedback *feedback )
271 : : {
272 : : // Check if the layer/mesh is valid
273 : 0 : if ( !mMeshLayer )
274 : 0 : return QgsGeometry();
275 : :
276 : 0 : populateCache( index, method );
277 : :
278 : 0 : return exportLines( value, feedback );
279 : 0 : }
280 : :
281 : 0 : QgsGeometry QgsMeshContours::exportLines( double value, QgsFeedback *feedback )
282 : : {
283 : 0 : std::unique_ptr<QgsMultiLineString> multiLineString( new QgsMultiLineString() );
284 : 0 : QSet<QPair<int, int>> exactEdges;
285 : :
286 : : // STEP 1: Get Data
287 : 0 : QVector<QgsMeshVertex> vertices = mTriangularMesh.vertices();
288 : 0 : const QVector<int> &trianglesToNativeFaces = mTriangularMesh.trianglesToNativeFaces();
289 : :
290 : : // STEP 2: For each triangle get the contour line
291 : 0 : for ( int i = 0; i < mTriangularMesh.triangles().size(); ++i )
292 : : {
293 : 0 : if ( feedback && feedback->isCanceled() )
294 : 0 : break;
295 : :
296 : 0 : int nativeIndex = trianglesToNativeFaces.at( i );
297 : 0 : if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
298 : 0 : continue;
299 : :
300 : 0 : const QgsMeshFace &triangle = mTriangularMesh.triangles().at( i );
301 : :
302 : : const int indices[3] =
303 : 0 : {
304 : 0 : triangle.at( 0 ),
305 : 0 : triangle.at( 1 ),
306 : 0 : triangle.at( 2 )
307 : : };
308 : :
309 : : const QVector<QgsMeshVertex> coords =
310 : 0 : {
311 : 0 : vertices.at( indices[0] ),
312 : 0 : vertices.at( indices[1] ),
313 : 0 : vertices.at( indices[2] )
314 : : };
315 : :
316 : : const double values[3] =
317 : 0 : {
318 : 0 : mDatasetValues.at( indices[0] ),
319 : 0 : mDatasetValues.at( indices[1] ),
320 : 0 : mDatasetValues.at( indices[2] )
321 : : };
322 : :
323 : : // any value is NaN
324 : 0 : if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
325 : 0 : continue;
326 : :
327 : : // value is outside the range
328 : 0 : if ( ( ( value > values[0] ) && ( value > values[1] ) && ( value > values[2] ) ) ||
329 : 0 : ( ( value < values[0] ) && ( value < values[1] ) && ( value < values[2] ) ) )
330 : 0 : continue;
331 : :
332 : : // all values are the same
333 : 0 : if ( qgsDoubleNear( values[0], values[1] ) && qgsDoubleNear( values[1], values[2] ) )
334 : 0 : continue;
335 : :
336 : : // go through all edges
337 : 0 : QgsPoint tmp;
338 : :
339 : 0 : for ( int i = 0; i < 3; ++i )
340 : : {
341 : 0 : const int j = ( i + 1 ) % 3;
342 : : // value is outside the range
343 : 0 : if ( ( ( value > values[i] ) && ( value > values[j] ) ) ||
344 : 0 : ( ( value < values[i] ) && ( value < values[j] ) ) )
345 : 0 : continue;
346 : :
347 : : // the whole edge is result and we are done
348 : 0 : if ( qgsDoubleNear( values[i], values[j] ) && qgsDoubleNear( values[i], values[j] ) )
349 : : {
350 : 0 : if ( exactEdges.contains( { indices[i], indices[j] } ) || exactEdges.contains( { indices[j], indices[i] } ) )
351 : : {
352 : 0 : break;
353 : : }
354 : : else
355 : : {
356 : 0 : exactEdges.insert( { indices[i], indices[j] } );
357 : 0 : std::unique_ptr<QgsLineString> line( new QgsLineString( coords[i], coords[j] ) );
358 : 0 : multiLineString->addGeometry( line.release() );
359 : : break;
360 : 0 : }
361 : : }
362 : :
363 : : // only one point matches, we are not interested in this
364 : 0 : if ( qgsDoubleNear( values[i], value ) || qgsDoubleNear( values[j], value ) )
365 : : {
366 : 0 : continue;
367 : : }
368 : :
369 : : // ok part of the result contour line is one point on this edge
370 : 0 : const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
371 : 0 : const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
372 : :
373 : 0 : if ( std::isnan( tmp.x() ) )
374 : : {
375 : : // ok we have found start point of the contour line
376 : 0 : tmp = xy;
377 : 0 : }
378 : : else
379 : : {
380 : : // we have found the end point of the contour line, we are done
381 : 0 : std::unique_ptr<QgsLineString> line( new QgsLineString( tmp, xy ) );
382 : 0 : multiLineString->addGeometry( line.release() );
383 : : break;
384 : 0 : }
385 : 0 : }
386 : 0 : }
387 : :
388 : : // STEP 3: merge the contour segments to linestrings
389 : 0 : if ( multiLineString->isEmpty() )
390 : : {
391 : 0 : return QgsGeometry();
392 : : }
393 : : else
394 : : {
395 : 0 : const QgsGeometry in( multiLineString.release() );
396 : 0 : const QgsGeometry res = in.mergeLines();
397 : 0 : return res;
398 : 0 : }
399 : 0 : }
400 : :
401 : 0 : void QgsMeshContours::populateCache( const QgsMeshDatasetIndex &index, QgsMeshRendererScalarSettings::DataResamplingMethod method )
402 : : {
403 : 0 : if ( !mMeshLayer )
404 : 0 : return;
405 : :
406 : 0 : if ( mCachedIndex != index )
407 : : {
408 : 0 : bool scalarDataOnVertices = mMeshLayer->dataProvider()->datasetGroupMetadata( index ).dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
409 : 0 : int count = scalarDataOnVertices ? mNativeMesh.vertices.count() : mNativeMesh.faces.count();
410 : :
411 : : // populate scalar values
412 : 0 : QgsMeshDataBlock vals = QgsMeshLayerUtils::datasetValues(
413 : 0 : mMeshLayer,
414 : 0 : index,
415 : : 0,
416 : 0 : count );
417 : 0 : if ( vals.isValid() )
418 : : {
419 : : // vals could be scalar or vectors, for contour rendering we want always magnitude
420 : 0 : mDatasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
421 : 0 : }
422 : : else
423 : : {
424 : 0 : mDatasetValues = QVector<double>( count, std::numeric_limits<double>::quiet_NaN() );
425 : : }
426 : :
427 : : // populate face active flag, always defined on faces
428 : 0 : mScalarActiveFaceFlagValues = mMeshLayer->dataProvider()->areFacesActive(
429 : 0 : index,
430 : : 0,
431 : 0 : mNativeMesh.faces.count() );
432 : :
433 : : // for data on faces, there could be request to interpolate the data to vertices
434 : 0 : if ( ( !scalarDataOnVertices ) )
435 : : {
436 : 0 : mDatasetValues = QgsMeshLayerUtils::interpolateFromFacesData(
437 : 0 : mDatasetValues,
438 : 0 : &mNativeMesh,
439 : 0 : &mTriangularMesh,
440 : 0 : &mScalarActiveFaceFlagValues,
441 : 0 : method
442 : : );
443 : 0 : }
444 : 0 : }
445 : 0 : }
|
