QGIS API Documentation 3.39.0-Master (be2050b798e)
Searching...
No Matches
qgstriangularmesh.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgstriangularmesh.cpp
3 ---------------------
4 begin : April 2018
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 *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
18#include <memory>
19#include <QList>
20#include "qgspolygon.h"
21#include "qgslinestring.h"
22#include "qgstriangularmesh.h"
23#include "qgsrendercontext.h"
25#include "qgsgeometry.h"
26#include "qgsrectangle.h"
27#include "qgslogger.h"
28#include "qgsmeshspatialindex.h"
29#include "qgsmeshlayerutils.h"
30#include "meshOptimizer/meshoptimizer.h"
31
32static void ENP_centroid_step( const QPolygonF &pX, double &cx, double &cy, double &signedArea, int i, int i1 )
33{
34 double x0 = 0.0; // Current vertex X
35 double y0 = 0.0; // Current vertex Y
36 double x1 = 0.0; // Next vertex X
37 double y1 = 0.0; // Next vertex Y
38 double a = 0.0; // Partial signed area
39
40 x0 = pX[i].x();
41 y0 = pX[i].y();
42 x1 = pX[i1].x();
43 y1 = pX[i1].y();
44 a = x0 * y1 - x1 * y0;
45 signedArea += a;
46 cx += ( x0 + x1 ) * a;
47 cy += ( y0 + y1 ) * a;
48}
49
50static void ENP_centroid( const QPolygonF &pX, double &cx, double &cy )
51{
52 // http://stackoverflow.com/questions/2792443/finding-the-centroid-of-a-polygon/2792459#2792459
53 cx = 0;
54 cy = 0;
55
56 if ( pX.isEmpty() )
57 return;
58
59 double signedArea = 0.0;
60
61 const QPointF &pt0 = pX.first();
62 QPolygonF localPolygon( pX.count() );
63 for ( int i = 0; i < pX.count(); ++i )
64 localPolygon[i] = pX.at( i ) - pt0;
65
66 // For all vertices except last
67 int i = 0;
68 for ( ; i < localPolygon.size() - 1; ++i )
69 {
70 ENP_centroid_step( localPolygon, cx, cy, signedArea, i, i + 1 );
71 }
72 // Do last vertex separately to avoid performing an expensive
73 // modulus operation in each iteration.
74 ENP_centroid_step( localPolygon, cx, cy, signedArea, i, 0 );
75
76 signedArea *= 0.5;
77 cx /= ( 6.0 * signedArea );
78 cy /= ( 6.0 * signedArea );
79
80 cx = cx + pt0.x();
81 cy = cy + pt0.y();
82}
83
84static void triangulateFaces( const QgsMeshFace &face,
85 int nativeIndex,
86 QVector<QgsMeshFace> &destinationFaces,
87 QVector<int> &triangularToNative,
88 const QgsMesh &verticesMeshSource )
89{
90 int vertexCount = face.size();
91 if ( vertexCount < 3 )
92 return;
93
94 while ( vertexCount > 3 )
95 {
96 // clip one ear from last 2 and first vertex
97 const QgsMeshFace ear = { face[vertexCount - 2], face[vertexCount - 1], face[0] };
98 if ( !( std::isnan( verticesMeshSource.vertex( ear[0] ).x() ) ||
99 std::isnan( verticesMeshSource.vertex( ear[1] ).x() ) ||
100 std::isnan( verticesMeshSource.vertex( ear[2] ).x() ) ) )
101 {
102 destinationFaces.push_back( ear );
103 triangularToNative.push_back( nativeIndex );
104 }
105 --vertexCount;
106 }
107
108 const QgsMeshFace triangle = { face[1], face[2], face[0] };
109 if ( !( std::isnan( verticesMeshSource.vertex( triangle[0] ).x() ) ||
110 std::isnan( verticesMeshSource.vertex( triangle[1] ).x() ) ||
111 std::isnan( verticesMeshSource.vertex( triangle[2] ).x() ) ) )
112 {
113 destinationFaces.push_back( triangle );
114 triangularToNative.push_back( nativeIndex );
115 }
116}
117
118void QgsTriangularMesh::triangulate( const QgsMeshFace &face, int nativeIndex )
119{
120 triangulateFaces( face, nativeIndex, mTriangularMesh.faces, mTrianglesToNativeFaces, mTriangularMesh );
121}
122
123QgsMeshVertex QgsTriangularMesh::transformVertex( const QgsMeshVertex &vertex, Qgis::TransformDirection direction ) const
124{
125 QgsMeshVertex transformedVertex = vertex;
126
127 if ( mCoordinateTransform.isValid() )
128 {
129 try
130 {
131 transformedVertex.transform( mCoordinateTransform, direction );
132 }
133 catch ( QgsCsException &cse )
134 {
135 Q_UNUSED( cse )
136 QgsDebugError( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ) );
137 transformedVertex = QgsMeshVertex();
138 }
139 }
140
141 return transformedVertex;
142}
143
144QgsMeshVertex QgsTriangularMesh::calculateCentroid( const QgsMeshFace &nativeFace )
145{
146 return QgsMeshUtils::centroid( nativeFace, mTriangularMesh.vertices );
147}
148
150{
151 return mAverageTriangleSize;
152}
153
156
157bool QgsTriangularMesh::update( QgsMesh *nativeMesh, const QgsCoordinateTransform &transform )
158{
159 Q_ASSERT( nativeMesh );
160
161 bool needUpdateVerticesCoordinates = mTriangularMesh.vertices.size() != nativeMesh->vertices.size() ||
162 ( ( mCoordinateTransform.isValid() || transform.isValid() ) &&
163 ( mCoordinateTransform.sourceCrs() != transform.sourceCrs() ||
164 mCoordinateTransform.destinationCrs() != transform.destinationCrs() ||
165 mCoordinateTransform.isValid() != transform.isValid() ) ) ;
166
167 bool needUpdateFrame = mTriangularMesh.vertices.size() != nativeMesh->vertices.size() ||
168 mNativeMeshFaceCentroids.size() != nativeMesh->faces.size() ||
169 mTriangularMesh.faces.size() < nativeMesh->faces.size() ||
170 mTriangularMesh.edges.size() != nativeMesh->edges.size();
171
172
173 // FIND OUT IF UPDATE IS NEEDED
174 if ( ! needUpdateVerticesCoordinates && !needUpdateFrame )
175 return false;
176
177 // CLEAN-UP
178 mTriangularMesh.vertices.clear();
179 if ( needUpdateFrame )
180 {
181 mTriangularMesh.faces.clear();
182 mTriangularMesh.edges.clear();
183 mEdgesToNativeEdges.clear();
184 mTrianglesToNativeFaces.clear();
185 }
186
187 // TRANSFORM VERTICES
188 mCoordinateTransform = transform;
189 mTriangularMesh.vertices.resize( nativeMesh->vertices.size() );
190 mExtent.setNull();
191 for ( int i = 0; i < nativeMesh->vertices.size(); ++i )
192 {
193 mTriangularMesh.vertices[i] = nativeToTriangularCoordinates( nativeMesh->vertices.at( i ) );
194 mExtent.include( mTriangularMesh.vertices.at( i ) );
195 }
196
197 if ( needUpdateFrame )
198 {
199 // CREATE TRIANGULAR MESH
200 for ( int i = 0; i < nativeMesh->faces.size(); ++i )
201 {
202 const QgsMeshFace &face = nativeMesh->faces.at( i ) ;
203 triangulate( face, i );
204 }
205 }
206
207 // CALCULATE CENTROIDS
208 mNativeMeshFaceCentroids.resize( nativeMesh->faces.size() );
209 for ( int i = 0; i < nativeMesh->faces.size(); ++i )
210 {
211 const QgsMeshFace &face = nativeMesh->faces.at( i ) ;
212 mNativeMeshFaceCentroids[i] = calculateCentroid( face );
213 }
214
215 // CALCULATE SPATIAL INDEX
216 mSpatialFaceIndex = QgsMeshSpatialIndex( mTriangularMesh, nullptr, QgsMesh::ElementType::Face );
217
218 if ( needUpdateFrame )
219 {
220 // SET ALL TRIANGLE CCW AND COMPUTE AVERAGE SIZE
221 finalizeTriangles();
222 }
223
224 // CREATE EDGES
225 // remove all edges with invalid vertices
226 if ( needUpdateFrame )
227 {
228 const QVector<QgsMeshEdge> edges = nativeMesh->edges;
229 for ( int nativeIndex = 0; nativeIndex < edges.size(); ++nativeIndex )
230 {
231 const QgsMeshEdge &edge = edges.at( nativeIndex );
232 if ( !( std::isnan( mTriangularMesh.vertex( edge.first ).x() ) ||
233 std::isnan( mTriangularMesh.vertex( edge.second ).x() ) ) )
234 {
235 mTriangularMesh.edges.push_back( edge );
236 mEdgesToNativeEdges.push_back( nativeIndex );
237 }
238 }
239 }
240
241 // CALCULATE CENTROIDS
242 mNativeMeshEdgeCentroids.resize( nativeMesh->edgeCount() );
243 for ( int i = 0; i < nativeMesh->edgeCount(); ++i )
244 {
245 const QgsMeshEdge &edge = nativeMesh->edges.at( i ) ;
246 const QgsPoint &a = mTriangularMesh.vertices[edge.first];
247 const QgsPoint &b = mTriangularMesh.vertices[edge.second];
248 mNativeMeshEdgeCentroids[i] = QgsMeshVertex( ( a.x() + b.x() ) / 2.0, ( a.y() + b.y() ) / 2.0, ( a.z() + b.z() ) / 2.0 );
249 }
250
251 // CALCULATE SPATIAL INDEX
252 mSpatialEdgeIndex = QgsMeshSpatialIndex( mTriangularMesh, nullptr, QgsMesh::ElementType::Edge );
253
254 return true;
255}
256
258{
259 return update( nativeMesh, mCoordinateTransform );
260}
261
262void QgsTriangularMesh::finalizeTriangles()
263{
264 mAverageTriangleSize = 0;
265 for ( int i = 0; i < mTriangularMesh.faceCount(); ++i )
266 {
267 QgsMeshFace &face = mTriangularMesh.faces[i];
268
269 const QgsMeshVertex &v0 = mTriangularMesh.vertex( face[0] );
270 const QgsMeshVertex &v1 = mTriangularMesh.vertex( face[1] );
271 const QgsMeshVertex &v2 = mTriangularMesh.vertex( face[2] );
272
273 QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox( v0, v1, v2 );
274
275 mAverageTriangleSize += std::fmax( bbox.width(), bbox.height() );
276
277 QgsMeshUtils::setCounterClockwise( face, v0, v1, v2 );
278 }
279 mAverageTriangleSize /= mTriangularMesh.faceCount();
280}
281
283{
284 return transformVertex( vertex, Qgis::TransformDirection::Forward );
285}
286
288{
289 return transformVertex( vertex, Qgis::TransformDirection::Reverse );
290}
291
293{
295 if ( !mCoordinateTransform.isShortCircuited() )
296 {
297 try
298 {
299 QgsCoordinateTransform extentTransform = mCoordinateTransform;
300 extentTransform.setBallparkTransformsAreAppropriate( true );
302 }
303 catch ( QgsCsException &cse )
304 {
305 Q_UNUSED( cse )
306 QgsDebugError( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ) );
307 }
308 }
309 else
311
312 return nativeExtent;
313}
314
316{
317 if ( !mIsExtentValid )
318 {
319 mExtent.setNull();
320 for ( int i = 0; i < mTriangularMesh.vertices.size(); ++i )
321 if ( !mTriangularMesh.vertices.at( i ).isEmpty() )
322 mExtent.include( mTriangularMesh.vertices.at( i ) );
323
324 mIsExtentValid = true;
325 }
326 return mExtent;
327}
328
330{
331 return mLod;
332}
333
335{
336 switch ( type )
337 {
339 return mTriangularMesh.vertexCount() != 0;
341 return mTriangularMesh.edgeCount() != 0;
343 return mTriangularMesh.faceCount() != 0;
344 }
345
346 return false;
347}
348
349void QgsTriangularMesh::addVertex( const QgsMeshVertex &vertex )
350{
351 QgsMeshVertex vertexInTriangularCoordinates = nativeToTriangularCoordinates( vertex );
352 mTriangularMesh.vertices.append( vertexInTriangularCoordinates );
353 if ( !vertexInTriangularCoordinates.isEmpty() )
354 mExtent.include( vertexInTriangularCoordinates );
355}
356
357const QVector<QgsMeshVertex> &QgsTriangularMesh::vertices() const
358{
359 return mTriangularMesh.vertices;
360}
361
362const QVector<QgsMeshFace> &QgsTriangularMesh::triangles() const
363{
364 return mTriangularMesh.faces;
365}
366
367const QVector<QgsMeshEdge> &QgsTriangularMesh::edges() const
368{
369 return mTriangularMesh.edges;
370}
371
372const QVector<QgsMeshVertex> &QgsTriangularMesh::centroids() const
373{
374 return faceCentroids();
375}
376
377const QVector<QgsMeshVertex> &QgsTriangularMesh::faceCentroids() const
378{
379 return mNativeMeshFaceCentroids;
380}
381
382const QVector<QgsMeshVertex> &QgsTriangularMesh::edgeCentroids() const
383{
384 return mNativeMeshEdgeCentroids;
385}
386
388{
389 return mTrianglesToNativeFaces;
390}
391
392const QVector<int> &QgsTriangularMesh::edgesToNativeEdges() const
393{
394 return mEdgesToNativeEdges;
395}
396
398{
399 QgsPointXY mapPoint;
400 if ( mCoordinateTransform.isValid() )
401 {
402 try
403 {
404 mapPoint = mCoordinateTransform.transform( point );
405 }
406 catch ( QgsCsException &cse )
407 {
408 Q_UNUSED( cse )
409 QgsDebugError( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ) );
410 mapPoint = point;
411 }
412 }
413 else
414 mapPoint = point;
415
416 return point;
417}
418
420{
421 const QList<int> faceIndexes = mSpatialFaceIndex.intersects( QgsRectangle( point, point ) );
422 for ( const int faceIndex : faceIndexes )
423 {
424 const QgsMeshFace &face = mTriangularMesh.faces.at( faceIndex );
425 const QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
426 if ( geom.contains( &point ) )
427 return faceIndex;
428 }
429 return -1;
430}
431
433{
434 int triangleIndex = faceIndexForPoint_v2( point );
435 if ( triangleIndex == -1 )
436 return -1;
437
438 if ( triangleIndex < mTrianglesToNativeFaces.count() )
439 return mTrianglesToNativeFaces.at( triangleIndex );
440
441 return -1;
442}
443
445{
446 QSet<int> concernedFaceIndex = QgsMeshUtils::nativeFacesFromTriangles(
447 faceIndexesForRectangle( rectangle ),
449 return concernedFaceIndex.values();
450}
451
453{
454 const QList<int> faceIndexes = mSpatialFaceIndex.intersects( QgsRectangle( point, point ) );
455
456 for ( const int faceIndex : faceIndexes )
457 {
458 const QgsMeshFace &face = mTriangularMesh.faces.at( faceIndex );
459 if ( QgsMeshUtils::isInTriangleFace( point, face, mTriangularMesh.vertices ) )
460 return faceIndex;
461 }
462 return -1;
463}
464
466{
467 return mSpatialFaceIndex.intersects( rectangle );
468}
469
471{
472 return mSpatialEdgeIndex.intersects( rectangle );
473}
474
475QVector<QVector3D> QgsTriangularMesh::vertexNormals( float vertScale ) const
476{
477 QVector<QVector3D> normals( vertices().count(), QVector3D( 0, 0, 0 ) );
478
479 for ( const auto &face : triangles() )
480 {
481 if ( face.isEmpty() )
482 continue;
483
484 for ( int i = 0; i < 3; i++ )
485 {
486 int index1( face.at( i ) );
487 int index2( face.at( ( i + 1 ) % 3 ) );
488 int index3( face.at( ( i + 2 ) % 3 ) );
489
490 const QgsMeshVertex &vert( vertices().at( index1 ) );
491 const QgsMeshVertex &otherVert1( vertices().at( index2 ) );
492 const QgsMeshVertex &otherVert2( vertices().at( index3 ) );
493
494 QVector3D v1( float( otherVert1.x() - vert.x() ), float( otherVert1.y() - vert.y() ), vertScale * float( otherVert1.z() - vert.z() ) );
495 QVector3D v2( float( otherVert2.x() - vert.x() ), float( otherVert2.y() - vert.y() ), vertScale * float( otherVert2.z() - vert.z() ) );
496
497 normals[index1] += QVector3D::crossProduct( v1, v2 );
498 }
499 }
500 return normals;
501}
502
503QVector<QgsTriangularMesh *> QgsTriangularMesh::simplifyMesh( double reductionFactor, int minimumTrianglesCount ) const
504{
505 QVector<QgsTriangularMesh *> simplifiedMeshes;
506
507 if ( mTriangularMesh.edgeCount() != 0 )
508 return simplifiedMeshes;
509
510 if ( !( reductionFactor > 1 ) )
511 return simplifiedMeshes;
512
513 size_t verticesCount = size_t( mTriangularMesh.vertices.count() );
514
515 unsigned int baseIndexCount = mTriangularMesh.faceCount() * 3;
516
517 QVector<unsigned int> indexes( mTriangularMesh.faces.count() * 3 );
518 for ( int i = 0; i < mTriangularMesh.faceCount(); ++i )
519 {
520 const QgsMeshFace &f = mTriangularMesh.face( i );
521 for ( int j = 0; j < 3; ++j )
522 indexes[i * 3 + j] = f.at( j );
523 }
524
525 QVector<float> vertices( mTriangularMesh.vertices.count() * 3 );
526 for ( int i = 0; i < mTriangularMesh.vertices.count(); ++i )
527 {
528 const QgsMeshVertex &v = mTriangularMesh.vertex( i );
529 vertices[i * 3] = v.x() ;
530 vertices[i * 3 + 1] = v.y() ;
531 vertices[i * 3 + 2] = v.z() ;
532 }
533
534 int path = 0;
535 while ( true )
536 {
537 QgsTriangularMesh *simplifiedMesh = new QgsTriangularMesh( *this );
538 size_t maxNumberOfIndexes = baseIndexCount / pow( reductionFactor, path + 1 );
539
540 if ( indexes.size() <= int( maxNumberOfIndexes ) )
541 {
542 delete simplifiedMesh;
543 break;
544 }
545
546 QVector<unsigned int> returnIndexes( indexes.size() );
547 //returned size could be different than goal size but not than the input indexes count
548 size_t size = meshopt_simplifySloppy(
549 returnIndexes.data(),
550 indexes.data(),
551 indexes.size(),
552 vertices.data(),
553 verticesCount,
554 sizeof( float ) * 3,
555 maxNumberOfIndexes );
556
557
558 returnIndexes.resize( size );
559
560 if ( size == 0 || int( size ) >= indexes.size() )
561 {
562 QgsDebugError( QStringLiteral( "Mesh simplification failed after %1 path" ).arg( path + 1 ) );
563 delete simplifiedMesh;
564 break;
565 }
566
567 QgsMesh newMesh;
568 newMesh.vertices = mTriangularMesh.vertices;
569
570 newMesh.faces.resize( returnIndexes.size() / 3 );
571 for ( int i = 0; i < newMesh.faces.size(); ++i )
572 {
573 QgsMeshFace f( 3 );
574 for ( size_t j = 0; j < 3 ; ++j )
575 f[j] = returnIndexes.at( i * 3 + j ) ;
576 newMesh.faces[i ] = f;
577 }
578
579 simplifiedMesh->mTriangularMesh = newMesh;
580 simplifiedMesh->mSpatialFaceIndex = QgsMeshSpatialIndex( simplifiedMesh->mTriangularMesh );
581 simplifiedMesh->finalizeTriangles();
582 simplifiedMeshes.push_back( simplifiedMesh );
583
584 QgsDebugMsgLevel( QStringLiteral( "Simplified mesh created with %1 triangles" ).arg( newMesh.faceCount() ), 2 );
585
586 simplifiedMesh->mTrianglesToNativeFaces = QVector<int>( simplifiedMesh->triangles().count(), 0 );
587 for ( int i = 0; i < simplifiedMesh->mTrianglesToNativeFaces.count(); ++i )
588 {
589 QgsMeshFace triangle = simplifiedMesh->triangles().at( i );
590 double x = 0;
591 double y = 0;
592 for ( size_t j = 0; j < 3 ; ++j )
593 {
594 x += mTriangularMesh.vertex( triangle[j] ).x();
595 y += mTriangularMesh.vertex( triangle[j] ).y();
596 }
597 x /= 3;
598 y /= 3;
599 QgsPoint centroid( x, y );
600 int indexInBaseMesh = faceIndexForPoint_v2( centroid );
601
602 if ( indexInBaseMesh == -1 )
603 {
604 // sometime the centroid of simplified mesh could be outside the base mesh,
605 // so try with vertices of the simplified triangle
606 int j = 0;
607 while ( indexInBaseMesh == -1 && j < 3 )
608 indexInBaseMesh = faceIndexForPoint_v2( mTriangularMesh.vertex( triangle[j++] ) );
609 }
610
611 if ( indexInBaseMesh > -1 && indexInBaseMesh < mTrianglesToNativeFaces.count() )
612 simplifiedMesh->mTrianglesToNativeFaces[i] = mTrianglesToNativeFaces[indexInBaseMesh];
613 }
614
615 simplifiedMesh->mLod = path + 1;
616 simplifiedMesh->mBaseTriangularMesh = this;
617
618 if ( simplifiedMesh->triangles().count() < minimumTrianglesCount )
619 break;
620
621 indexes = returnIndexes;
622 ++path;
623 }
624
625 return simplifiedMeshes;
626}
627
628std::unique_ptr< QgsPolygon > QgsMeshUtils::toPolygon( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
629{
630 QVector<QgsPoint> ring;
631 for ( int j = 0; j < face.size(); ++j )
632 {
633 int vertexId = face[j];
634 Q_ASSERT( vertexId < vertices.size() );
635 const QgsPoint &vertex = vertices[vertexId];
636 ring.append( vertex );
637 }
638 std::unique_ptr< QgsPolygon > polygon = std::make_unique< QgsPolygon >();
639 polygon->setExteriorRing( new QgsLineString( ring ) );
640 return polygon;
641}
642
643QgsGeometry QgsMeshUtils::toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
644{
645 return QgsGeometry( QgsMeshUtils::toPolygon( face, vertices ) );
646}
647
648static QSet<int> _nativeElementsFromElements( const QList<int> &indexes, const QVector<int> &elementToNativeElements )
649{
650 QSet<int> nativeElements;
651 for ( const int index : indexes )
652 {
653 if ( index < elementToNativeElements.count() )
654 {
655 const int nativeIndex = elementToNativeElements[index];
656 nativeElements.insert( nativeIndex );
657 }
658 }
659 return nativeElements;
660}
661
662QSet<int> QgsMeshUtils::nativeFacesFromTriangles( const QList<int> &triangleIndexes, const QVector<int> &trianglesToNativeFaces )
663{
664 return _nativeElementsFromElements( triangleIndexes, trianglesToNativeFaces );
665}
666
667QSet<int> QgsMeshUtils::nativeEdgesFromEdges( const QList<int> &edgesIndexes, const QVector<int> &edgesToNativeEdges )
668{
669 return _nativeElementsFromElements( edgesIndexes, edgesToNativeEdges );
670}
671
672
673static double _isLeft2D( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p )
674{
675 return ( p2.x() - p1.x() ) * ( p.y() - p1.y() ) - ( p.x() - p1.x() ) * ( p2.y() - p1.y() );
676}
677
678static bool _isInTriangle2D( const QgsPoint &p, const QVector<QgsMeshVertex> &triangle )
679{
680 return ( ( _isLeft2D( triangle[2], triangle[0], p ) * _isLeft2D( triangle[2], triangle[0], triangle[1] ) >= 0 )
681 && ( _isLeft2D( triangle[0], triangle[1], p ) * _isLeft2D( triangle[0], triangle[1], triangle[2] ) >= 0 )
682 && ( _isLeft2D( triangle[2], triangle[1], p ) * _isLeft2D( triangle[2], triangle[1], triangle[0] ) >= 0 ) );
683}
684
685bool QgsMeshUtils::isInTriangleFace( const QgsPointXY point, const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
686{
687 if ( face.count() != 3 )
688 return false;
689
690 QVector<QgsMeshVertex> triangle( 3 );
691 for ( int i = 0; i < 3; ++i )
692 {
693 if ( face[i] > vertices.count() )
694 return false;
695 triangle[i] = vertices[face[i]];
696 }
697
698 const QgsPoint p( point.x(), point.y() );
699
700 return _isInTriangle2D( p, triangle );
701}
702
703QSet<int> QgsMeshUtils::nativeVerticesFromTriangles( const QList<int> &triangleIndexes, const QVector<QgsMeshFace> &triangles )
704{
705 QSet<int> uniqueVertices;
706 for ( int triangleIndex : triangleIndexes )
707 {
708 const QgsMeshFace triangle = triangles[triangleIndex];
709 for ( int i : triangle )
710 {
711 uniqueVertices.insert( i );
712 }
713 }
714 return uniqueVertices;
715}
716
717QSet<int> QgsMeshUtils::nativeVerticesFromEdges( const QList<int> &edgesIndexes, const QVector<QgsMeshEdge> &edges )
718{
719 QSet<int> uniqueVertices;
720 for ( int edgeIndex : edgesIndexes )
721 {
722 const QgsMeshEdge edge = edges[edgeIndex];
723 uniqueVertices.insert( edge.first );
724 uniqueVertices.insert( edge.second );
725 }
726 return uniqueVertices;
727}
728
730{
731 //if necessary defined removes triangles index
732 if ( changes.mRemovedTriangleIndexes.isEmpty() && !changes.mNativeFaceIndexesToRemove.isEmpty() )
733 {
734 for ( int nf = 0; nf < changes.mNativeFaceIndexesToRemove.count(); ++nf )
735 {
736 int nativeIndex = changes.mNativeFaceIndexesToRemove.at( nf );
737 const QgsMeshFace &nativeFace = changes.mNativeFacesToRemove.at( nf );
738 Q_ASSERT( !nativeFace.isEmpty() );
739
740 QgsRectangle nativeFaceExtent( mTriangularMesh.vertex( nativeFace.at( 0 ) ), mTriangularMesh.vertex( nativeFace.at( 0 ) ) );
741 for ( int i = 1; i < nativeFace.count(); ++i )
742 {
743 const QgsMeshVertex &triangularVertex = mTriangularMesh.vertex( nativeFace.at( i ) );
744 nativeFaceExtent.include( triangularVertex );
745 }
746
747 QList<int> concernedTriangle = faceIndexesForRectangle( nativeFaceExtent );
748 //Remove only those corresponding to the native face
749 for ( int i = 0; i < concernedTriangle.count(); ++i )
750 {
751 int triangleIndex = concernedTriangle.at( i );
752 if ( mTrianglesToNativeFaces.at( triangleIndex ) == nativeIndex )
753 changes.mRemovedTriangleIndexes.append( triangleIndex );
754 }
755 }
756 }
757
758 if ( changes.mOldZValue.isEmpty() && !changes.mNewZValue.isEmpty() )
759 {
760 changes.mOldZValue.reserve( changes.mNewZValue.count() );
761 for ( int i = 0; i < changes.mNewZValue.count(); ++i )
762 changes.mOldZValue.append( mTriangularMesh.vertices.at( changes.mChangedVerticesCoordinates.at( i ) ).z() );
763 }
764
765 if ( changes.mTriangleIndexesGeometryChanged.isEmpty() && !changes.mNativeFaceIndexesGeometryChanged.isEmpty() )
766 {
767 for ( int i = 0; i < changes.mNativeFaceIndexesGeometryChanged.count(); ++i )
768 {
769 const QgsMeshFace &nativeFace = changes.mNativeFacesGeometryChanged.at( i );
770 if ( nativeFace.count() < 2 )
771 continue;
772 QgsRectangle bbox( mTriangularMesh.vertices.at( nativeFace.at( 0 ) ), mTriangularMesh.vertices.at( nativeFace.at( 1 ) ) );
773
774 for ( int i = 2; i < nativeFace.count(); ++i )
775 bbox.include( mTriangularMesh.vertices.at( nativeFace.at( i ) ) );
776
777 QList<int> triangleIndexes = faceIndexesForRectangle( bbox );
778 int pos = 0;
779 while ( pos < triangleIndexes.count() )
780 {
781 if ( trianglesToNativeFaces().at( triangleIndexes.at( pos ) ) !=
782 changes.mNativeFaceIndexesGeometryChanged.at( i ) )
783 triangleIndexes.removeAt( pos );
784 else
785 ++pos;
786 }
787 changes.mTriangleIndexesGeometryChanged.append( triangleIndexes );
788 }
789 }
790
792 for ( const QgsMeshVertex &vertex : std::as_const( changes.mAddedVertices ) )
794
797 {
799 int firstNewNativeFacesIndex = mNativeMeshFaceCentroids.count();
800 for ( int i = 0; i < changes.mNativeFacesToAdd.count(); ++i )
801 {
802 const QgsMeshFace &nativeFace = changes.mNativeFacesToAdd.at( i );
803 triangulate( nativeFace, firstNewNativeFacesIndex + i );
804 mNativeMeshFaceCentroids.append( calculateCentroid( nativeFace ) );
805 }
806
807 for ( int i = changes.mTrianglesAddedStartIndex; i < mTriangularMesh.faceCount(); ++i )
809 }
810
811 // remove faces
812 for ( int i = 0; i < changes.mRemovedTriangleIndexes.count(); ++i )
813 {
814 int triangleIndex = changes.mRemovedTriangleIndexes.at( i );
815 mTrianglesToNativeFaces[triangleIndex] = -1;
816 mSpatialFaceIndex.removeFace( triangleIndex, mTriangularMesh );
817 mTriangularMesh.faces[triangleIndex] = QgsMeshFace();
818 }
819
820 for ( int i = 0; i < changes.mNativeFaceIndexesToRemove.count(); ++i )
821 mNativeMeshFaceCentroids[changes.mNativeFaceIndexesToRemove.at( i )] = QgsMeshVertex();
822
823 // remove vertices
824 for ( int i = 0; i < changes.mVerticesIndexesToRemove.count(); ++i )
825 mTriangularMesh.vertices[changes.mVerticesIndexesToRemove.at( i )] = QgsMeshVertex();
826
827 if ( !changes.mVerticesIndexesToRemove.isEmpty() )
828 mIsExtentValid = false;
829
830 // change Z value
831 for ( int i = 0; i < changes.mNewZValue.count(); ++i )
832 {
833 int vertexIndex = changes.mChangedVerticesCoordinates.at( i );
834 mTriangularMesh.vertices[vertexIndex].setZ( changes.mNewZValue.at( i ) );
835 }
836
837 //remove outdated spatial index
838 for ( const int triangleIndex : std::as_const( changes.mTriangleIndexesGeometryChanged ) )
839 mSpatialFaceIndex.removeFace( triangleIndex, mTriangularMesh );
840
841 // change (X,Y) of vertices
842 for ( int i = 0; i < changes.mNewXYValue.count(); ++i )
843 {
844 const QgsPointXY &nativeCoordinates = changes.mNewXYValue.at( i );
845 const QgsMeshVertex nativeVertex( nativeCoordinates.x(),
846 nativeCoordinates.y(),
847 mTriangularMesh.vertices.at( changes.mChangedVerticesCoordinates.at( i ) ).z() );
848
849 mTriangularMesh.vertices[changes.mChangedVerticesCoordinates.at( i )] = nativeToTriangularCoordinates( nativeVertex );
850 }
851
852 //restore spatial undex
853 for ( const int triangleIndex : std::as_const( changes.mTriangleIndexesGeometryChanged ) )
855
856 //update native faces
857 for ( int i = 0; i < changes.mNativeFaceIndexesGeometryChanged.count(); ++i )
858 mNativeMeshFaceCentroids[changes.mNativeFaceIndexesGeometryChanged.at( i )] = calculateCentroid( changes.mNativeFacesGeometryChanged.at( i ) );
859}
860
862{
865 {
866 for ( int i = changes.mTrianglesAddedStartIndex; i < mTriangularMesh.faceCount(); ++i )
867 mSpatialFaceIndex.removeFace( i, mTriangularMesh );
868
869 int initialNativeFacesCount = mNativeMeshFaceCentroids.count() - changes.mNativeFacesToAdd.count();
870
873 mNativeMeshFaceCentroids.resize( initialNativeFacesCount );
874 }
875
876 int initialVerticesCount = mTriangularMesh.vertices.count() - changes.mAddedVertices.count();
877 mTriangularMesh.vertices.resize( initialVerticesCount );
878
880 mIsExtentValid = false;
881
882 // for each vertex to remove we need to update the vertices with the native vertex
883 for ( const int i : std::as_const( changes.mVerticesIndexesToRemove ) )
884 mTriangularMesh.vertices[i] = nativeToTriangularCoordinates( nativeMesh.vertex( i ) );
885
886 if ( !changes.mVerticesIndexesToRemove.isEmpty() )
887 mIsExtentValid = false;
888
889 // reverse removed faces
890 QVector<QgsMeshFace> restoredTriangles;
891 QVector<int> restoredTriangularToNative;
892 for ( int i = 0; i < changes.mNativeFacesToRemove.count(); ++i )
893 {
894 const QgsMeshFace &nativeFace = changes.mNativeFacesToRemove.at( i );
895 triangulateFaces( nativeFace,
896 changes.mNativeFaceIndexesToRemove.at( i ),
897 restoredTriangles,
898 restoredTriangularToNative,
899 mTriangularMesh );
900 mNativeMeshFaceCentroids[changes.mNativeFaceIndexesToRemove.at( i )] = calculateCentroid( nativeFace );
901 }
902 for ( int i = 0; i < changes.mRemovedTriangleIndexes.count(); ++i )
903 {
904 int triangleIndex = changes.mRemovedTriangleIndexes.at( i );
905 mTriangularMesh.faces[triangleIndex] = restoredTriangles.at( i );
907 mTrianglesToNativeFaces[triangleIndex] = restoredTriangularToNative.at( i );
908 }
909
910 // reverse Z value
911 for ( int i = 0; i < changes.mOldZValue.count(); ++i )
912 {
913 int vertexIndex = changes.mChangedVerticesCoordinates.at( i );
914 mTriangularMesh.vertices[vertexIndex].setZ( changes.mOldZValue.at( i ) );
915 }
916
917 //remove outdated spatial index
918 for ( const int triangleIndex : std::as_const( changes.mTriangleIndexesGeometryChanged ) )
919 mSpatialFaceIndex.removeFace( triangleIndex, mTriangularMesh );
920
921 // reverse (X,Y) of vertices
922 for ( int i = 0; i < changes.mOldXYValue.count(); ++i )
923 {
924 const QgsPointXY &nativeCoordinates = changes.mOldXYValue.at( i );
925 const QgsMeshVertex nativeVertex( nativeCoordinates.x(),
926 nativeCoordinates.y(),
927 mTriangularMesh.vertices.at( changes.mChangedVerticesCoordinates.at( i ) ).z() );
928
929 mTriangularMesh.vertices[changes.mChangedVerticesCoordinates.at( i )] = nativeToTriangularCoordinates( nativeVertex );
930 }
931
932 //restore spatial undex
933 for ( const int triangleIndex : std::as_const( changes.mTriangleIndexesGeometryChanged ) )
935
936 //update native faces
937 for ( int i = 0; i < changes.mNativeFaceIndexesGeometryChanged.count(); ++i )
938 mNativeMeshFaceCentroids[changes.mNativeFaceIndexesGeometryChanged.at( i )] = calculateCentroid( changes.mNativeFacesGeometryChanged.at( i ) );
939}
940
942 const QgsMesh &nativeMesh )
943{
945 mVerticesIndexesToRemove = topologicalChanges.verticesToRemoveIndexes();
947 mNativeFacesToRemove = topologicalChanges.removedFaces();
948 mNativeFaceIndexesToRemove = topologicalChanges.removedFaceIndexes();
949 mChangedVerticesCoordinates = topologicalChanges.changedCoordinatesVerticesIndexes();
950 mNewZValue = topologicalChanges.newVerticesZValues();
951 mNewXYValue = topologicalChanges.newVerticesXYValues();
952 mOldXYValue = topologicalChanges.oldVerticesXYValues();
953
954 mNativeFaceIndexesGeometryChanged = topologicalChanges.nativeFacesIndexesGeometryChanged();
955 mNativeFacesGeometryChanged.resize( mNativeFaceIndexesGeometryChanged.count() );
956 for ( int i = 0; i < mNativeFaceIndexesGeometryChanged.count(); ++i )
957 mNativeFacesGeometryChanged[i] = nativeMesh.face( mNativeFaceIndexesGeometryChanged.at( i ) );
958}
959
960QgsMeshVertex QgsMeshUtils::centroid( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
961{
962 QVector<QPointF> points( face.size() );
963 for ( int j = 0; j < face.size(); ++j )
964 {
965 int index = face.at( j );
966 const QgsMeshVertex &vertex = vertices.at( index ); // we need vertices in map coordinate
967 points[j] = vertex.toQPointF();
968 }
969 QPolygonF poly( points );
970 double cx, cy;
971 ENP_centroid( poly, cx, cy );
972 return QgsMeshVertex( cx, cy );
973}
974
975void QgsMeshUtils::setCounterClockwise( QgsMeshFace &triangle, const QgsMeshVertex &v0, const QgsMeshVertex &v1, const QgsMeshVertex &v2 )
976{
977 //To have consistent clock wise orientation of triangles which is necessary for 3D rendering
978 //Check the clock wise, and if it is not counter clock wise, swap indexes to make the oientation counter clock wise
979 double ux = v1.x() - v0.x();
980 double uy = v1.y() - v0.y();
981 double vx = v2.x() - v0.x();
982 double vy = v2.y() - v0.y();
983
984 double crossProduct = ux * vy - uy * vx;
985 if ( crossProduct < 0 ) //CW -->change the orientation
986 {
987 std::swap( triangle[1], triangle[2] );
988 }
989}
TransformDirection
Indicates the direction (forward or inverse) of a transform.
Definition qgis.h:2369
@ Forward
Forward transform (from source to destination)
@ Reverse
Reverse/inverse transform (from destination to source)
Class for doing transforms between two map coordinate systems.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from.
void setBallparkTransformsAreAppropriate(bool appropriate)
Sets whether approximate "ballpark" results are appropriate for this coordinate transform.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const
Transform the point from the source CRS to the destination CRS.
bool isShortCircuited() const
Returns true if the transform short circuits because the source and destination are equivalent.
QgsRectangle transformBoundingBox(const QgsRectangle &rectangle, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool handle180Crossover=false) const
Transforms a rectangle from the source CRS to the destination CRS.
bool isValid() const
Returns true if the coordinate transform is valid, ie both the source and destination CRS have been s...
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
Custom exception class for Coordinate Reference System related exceptions.
QString what() const
A geometry is the spatial representation of a feature.
bool contains(const QgsPointXY *p) const
Returns true if the geometry contains the point p.
Line string geometry type, with support for z-dimension and m-values.
A spatial index for QgsMeshFace or QgsMeshEdge objects.
QList< int > intersects(const QgsRectangle &rectangle) const
Returns a list of face ids with a bounding box which intersects the specified rectangle.
void addFace(int faceIndex, const QgsMesh &mesh)
Adds a face with faceIndex from the mesh in the spatial index.
void removeFace(int faceIndex, const QgsMesh &mesh)
Removes a face with faceIndex from the mesh in the spatial index.
A class to represent a 2D point.
Definition qgspointxy.h:60
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
QPointF toQPointF() const
Returns the point as a QPointF.
Definition qgspoint.h:382
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
bool isEmpty() const override
Returns true if the geometry is empty.
Definition qgspoint.cpp:737
void transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection d=Qgis::TransformDirection::Forward, bool transformZ=false) override
Transforms the geometry using a coordinate transform.
Definition qgspoint.cpp:383
double y
Definition qgspoint.h:53
A rectangle specified with double values.
void include(const QgsPointXY &p)
Updates the rectangle to include the specified point.
double width() const
Returns the width of the rectangle.
double height() const
Returns the height of the rectangle.
void setNull()
Mark a rectangle as being null (holding no spatial information).
Class that contains topological differences between two states of a topological mesh,...
QVector< QgsMeshFace > removedFaces() const
Returns the faces that are removed with this changes.
Returns the added vertices with this changes.
QList< int > changedCoordinatesVerticesIndexes() const
Returns the indexes of vertices that have changed coordinates.
QList< int > removedFaceIndexes() const
Returns the indexes of the faces that are removed with this changes.
QList< double > newVerticesZValues() const
Returns the new Z values of vertices that have changed their coordinates.
Returns the face that are added with this changes.
QList< QgsPointXY > oldVerticesXYValues() const
Returns the old (X,Y) values of vertices that have changed their coordinates.
QList< QgsPointXY > newVerticesXYValues() const
Returns the new (X,Y) values of vertices that have changed their coordinates.
QList< int > nativeFacesIndexesGeometryChanged() const
Returns a list of the native face indexes that have a geometry changed.
QList< int > verticesToRemoveIndexes() const
Returns the indexes of vertices to remove.
The Changes class is used to make changes of the triangular and to keep traces of this changes If a C...
Changes()=default
Default constructor, no changes.
Triangular/Derived Mesh is mesh with vertices in map coordinates.
const QVector< QgsMeshVertex > & edgeCentroids() const
Returns centroids of the native edges in map CRS.
const QVector< QgsMeshFace > & triangles() const
Returns triangles.
QVector< QgsTriangularMesh * > simplifyMesh(double reductionFactor, int minimumTrianglesCount=10) const
Returns simplified meshes.
QgsRectangle nativeExtent()
Returns the extent of the mesh in the native mesh coordinates system, returns empty extent if the tra...
QgsPointXY transformFromLayerToTrianglesCoordinates(const QgsPointXY &point) const
Transforms a point from layer coordinates system to triangular Mesh coordinates system.
QgsTriangularMesh()
Ctor.
int levelOfDetail() const
Returns the corresponding index of level of detail on which this mesh is associated.
QgsRectangle extent() const
Returns the extent of the triangular mesh in map coordinates.
int faceIndexForPoint(const QgsPointXY &point) const
Finds index of triangle at given point It uses spatial indexing.
int nativeFaceIndexForPoint(const QgsPointXY &point) const
Finds index of native face at given point It uses spatial indexing.
double averageTriangleSize() const
Returns the average size of triangles in map unit.
void reverseChanges(const Changes &changes, const QgsMesh &nativeMesh)
Reverses the changes on the triangular mesh (see Changes)
void applyChanges(const Changes &changes)
Applies the changes on the triangular mesh (see Changes)
QList< int > edgeIndexesForRectangle(const QgsRectangle &rectangle) const
Finds indexes of edges intersecting given bounding box It uses spatial indexing.
Q_DECL_DEPRECATED const QVector< QgsMeshVertex > & centroids() const
Returns centroids of the native faces in map CRS.
const QVector< QgsMeshVertex > & vertices() const
Returns vertices in map coordinate system.
const QVector< QgsMeshEdge > & edges() const
Returns edges.
bool contains(const QgsMesh::ElementType &type) const
Returns whether the mesh contains mesh elements of given type.
QgsMeshVertex triangularToNativeCoordinates(const QgsMeshVertex &vertex) const
Transforms the vertex from triangular mesh coordinates system to native coordinates system.
QVector< QVector3D > vertexNormals(float vertScale) const
Calculates and returns normal vector on each vertex that is part of any face.
QgsMeshVertex nativeToTriangularCoordinates(const QgsMeshVertex &vertex) const
Transforms the vertex from native coordinates system to triangular mesh coordinates system.
bool update(QgsMesh *nativeMesh, const QgsCoordinateTransform &transform)
Constructs triangular mesh from layer's native mesh and transform to destination CRS.
~QgsTriangularMesh()
Dtor.
const QVector< QgsMeshVertex > & faceCentroids() const
Returns centroids of the native faces in map CRS.
QList< int > nativeFaceIndexForRectangle(const QgsRectangle &rectangle) const
Finds indexes of native faces which bounding boxes intersect given bounding box It uses spatial index...
const QVector< int > & trianglesToNativeFaces() const
Returns mapping between triangles and original faces.
const QVector< int > & edgesToNativeEdges() const
Returns mapping between edges and original edges.
QList< int > faceIndexesForRectangle(const QgsRectangle &rectangle) const
Finds indexes of triangles intersecting given bounding box It uses spatial indexing.
int faceIndexForPoint_v2(const QgsPointXY &point) const
Finds index of triangle at given point It uses spatial indexing and don't use geos to be faster.
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
#define QgsDebugError(str)
Definition qgslogger.h:38
QVector< int > QgsMeshFace
List of vertex indexes.
QPair< int, int > QgsMeshEdge
Edge is a straight line seqment between 2 points.
QgsPoint QgsMeshVertex
xyz coords of vertex
Mesh - vertices, edges and faces.
int vertexCount() const
Returns number of vertices.
QVector< QgsMeshVertex > vertices
QgsMeshFace face(int index) const
Returns a face at the index.
QVector< QgsMeshFace > faces
int faceCount() const
Returns number of faces.
ElementType
Defines type of mesh elements.
QgsMeshVertex vertex(int index) const
Returns a vertex at the index.
int edgeCount() const
Returns number of edge.
QVector< QgsMeshEdge > edges