QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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 *
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 
18 #include <memory>
19 #include <QList>
20 #include "qgspolygon.h"
21 #include "qgslinestring.h"
22 #include "qgstriangularmesh.h"
23 #include "qgsrendercontext.h"
24 #include "qgscoordinatetransform.h"
25 #include "qgsgeometry.h"
26 #include "qgsrectangle.h"
27 #include "qgslogger.h"
28 #include "qgsmeshspatialindex.h"
29 
30 static void ENP_centroid_step( const QPolygonF &pX, double &cx, double &cy, double &signedArea, int i, int i1 )
31 {
32  double x0 = 0.0; // Current vertex X
33  double y0 = 0.0; // Current vertex Y
34  double x1 = 0.0; // Next vertex X
35  double y1 = 0.0; // Next vertex Y
36  double a = 0.0; // Partial signed area
37 
38  x0 = pX[i].x();
39  y0 = pX[i].y();
40  x1 = pX[i1].x();
41  y1 = pX[i1].y();
42  a = x0 * y1 - x1 * y0;
43  signedArea += a;
44  cx += ( x0 + x1 ) * a;
45  cy += ( y0 + y1 ) * a;
46 }
47 
48 static void ENP_centroid( const QPolygonF &pX, double &cx, double &cy )
49 {
50  // http://stackoverflow.com/questions/2792443/finding-the-centroid-of-a-polygon/2792459#2792459
51  cx = 0;
52  cy = 0;
53 
54  if ( pX.isEmpty() )
55  return;
56 
57  double signedArea = 0.0;
58 
59  // For all vertices except last
60  int i = 0;
61  for ( ; i < pX.size() - 1; ++i )
62  {
63  ENP_centroid_step( pX, cx, cy, signedArea, i, i + 1 );
64  }
65  // Do last vertex separately to avoid performing an expensive
66  // modulus operation in each iteration.
67  ENP_centroid_step( pX, cx, cy, signedArea, i, 0 );
68 
69  signedArea *= 0.5;
70  cx /= ( 6.0 * signedArea );
71  cy /= ( 6.0 * signedArea );
72 }
73 
74 
75 void QgsTriangularMesh::triangulate( const QgsMeshFace &face, int nativeIndex )
76 {
77  int vertexCount = face.size();
78  if ( vertexCount < 3 )
79  return;
80 
81  while ( vertexCount > 3 )
82  {
83  // clip one ear from last 2 and first vertex
84  const QgsMeshFace ear = { face[vertexCount - 2], face[vertexCount - 1], face[0] };
85  mTriangularMesh.faces.push_back( ear );
86  mTrianglesToNativeFaces.push_back( nativeIndex );
87  --vertexCount;
88  }
89 
90  const QgsMeshFace triangle = { face[1], face[2], face[0] };
91  mTriangularMesh.faces.push_back( triangle );
92  mTrianglesToNativeFaces.push_back( nativeIndex );
93 }
94 
97 
98 void QgsTriangularMesh::update( QgsMesh *nativeMesh, QgsRenderContext *context )
99 {
100  Q_ASSERT( nativeMesh );
101  Q_ASSERT( context );
102 
103  // FIND OUT IF UPDATE IS NEEDED
104  if ( mTriangularMesh.vertices.size() >= nativeMesh->vertices.size() &&
105  mTriangularMesh.faces.size() >= nativeMesh->faces.size() &&
106  mCoordinateTransform.sourceCrs() == context->coordinateTransform().sourceCrs() &&
107  mCoordinateTransform.destinationCrs() == context->coordinateTransform().destinationCrs() )
108  return;
109 
110  // CLEAN-UP
111  mTriangularMesh.vertices.clear();
112  mTriangularMesh.faces.clear();
113  mTrianglesToNativeFaces.clear();
114  mNativeMeshFaceCentroids.clear();
115 
116  // TRANSFORM VERTICES
117  mCoordinateTransform = context->coordinateTransform();
118  mTriangularMesh.vertices.resize( nativeMesh->vertices.size() );
119  for ( int i = 0; i < nativeMesh->vertices.size(); ++i )
120  {
121  const QgsMeshVertex &vertex = nativeMesh->vertices.at( i );
122  if ( mCoordinateTransform.isValid() )
123  {
124  try
125  {
126  QgsPointXY mapPoint = mCoordinateTransform.transform( QgsPointXY( vertex.x(), vertex.y() ) );
127  QgsMeshVertex mapVertex( mapPoint );
128  mapVertex.addZValue( vertex.z() );
129  mapVertex.setM( vertex.m() );
130  mTriangularMesh.vertices[i] = mapVertex;
131  }
132  catch ( QgsCsException &cse )
133  {
134  Q_UNUSED( cse )
135  QgsDebugMsg( QStringLiteral( "Caught CRS exception %1" ).arg( cse.what() ) );
136  mTriangularMesh.vertices[i] = vertex;
137  }
138  }
139  else
140  {
141  mTriangularMesh.vertices[i] = vertex;
142  }
143  }
144 
145  // CREATE TRIANGULAR MESH
146  for ( int i = 0; i < nativeMesh->faces.size(); ++i )
147  {
148  const QgsMeshFace &face = nativeMesh->faces.at( i ) ;
149  triangulate( face, i );
150  }
151 
152  // CALCULATE CENTROIDS
153  mNativeMeshFaceCentroids.resize( nativeMesh->faces.size() );
154  for ( int i = 0; i < nativeMesh->faces.size(); ++i )
155  {
156  const QgsMeshFace &face = nativeMesh->faces.at( i ) ;
157  QVector<QPointF> points;
158  points.reserve( face.size() );
159  for ( int j = 0; j < face.size(); ++j )
160  {
161  int index = face.at( j );
162  const QgsMeshVertex &vertex = mTriangularMesh.vertices[index]; // we need projected vertices
163  points.push_back( vertex.toQPointF() );
164  }
165  QPolygonF poly( points );
166  double cx, cy;
167  ENP_centroid( poly, cx, cy );
168  mNativeMeshFaceCentroids[i] = QgsMeshVertex( cx, cy );
169  }
170 
171  // CALCULATE SPATIAL INDEX
172  mSpatialIndex = QgsMeshSpatialIndex( mTriangularMesh );
173 }
174 
175 const QVector<QgsMeshVertex> &QgsTriangularMesh::vertices() const
176 {
177  return mTriangularMesh.vertices;
178 }
179 
180 const QVector<QgsMeshFace> &QgsTriangularMesh::triangles() const
181 {
182  return mTriangularMesh.faces;
183 }
184 
185 const QVector<QgsMeshVertex> &QgsTriangularMesh::centroids() const
186 {
187  return mNativeMeshFaceCentroids;
188 }
189 
190 const QVector<int> &QgsTriangularMesh::trianglesToNativeFaces() const
191 {
192  return mTrianglesToNativeFaces;
193 }
194 
196 {
197  const QList<int> faceIndexes = mSpatialIndex.intersects( QgsRectangle( point, point ) );
198  for ( const int faceIndex : faceIndexes )
199  {
200  const QgsMeshFace &face = mTriangularMesh.faces.at( faceIndex );
201  const QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
202  if ( geom.contains( &point ) )
203  return faceIndex;
204  }
205  return -1;
206 }
207 
208 QList<int> QgsTriangularMesh::faceIndexesForRectangle( const QgsRectangle &rectangle ) const
209 {
210  return mSpatialIndex.intersects( rectangle );
211 }
212 
213 std::unique_ptr< QgsPolygon > QgsMeshUtils::toPolygon( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
214 {
215  QVector<QgsPoint> ring;
216  for ( int j = 0; j < face.size(); ++j )
217  {
218  int vertexId = face[j];
219  Q_ASSERT( vertexId < vertices.size() );
220  const QgsPoint &vertex = vertices[vertexId];
221  ring.append( vertex );
222  }
223  std::unique_ptr< QgsPolygon > polygon = qgis::make_unique< QgsPolygon >();
224  polygon->setExteriorRing( new QgsLineString( ring ) );
225  return polygon;
226 }
227 
228 QgsGeometry QgsMeshUtils::toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
229 {
230  return QgsGeometry( QgsMeshUtils::toPolygon( face, vertices ) );
231 }
232 
233 QList<int> QgsMeshUtils::nativeFacesFromTriangles( const QList<int> &triangleIndexes, const QVector<int> &trianglesToNativeFaces )
234 {
235  QSet<int> nativeFaces;
236  for ( const int triangleIndex : triangleIndexes )
237  {
238  const int nativeIndex = trianglesToNativeFaces[triangleIndex];
239  nativeFaces.insert( nativeIndex );
240  }
241  return nativeFaces.toList();
242 }
CORE_EXPORT QgsGeometry toGeometry(const QgsMeshFace &face, const QVector< QgsMeshVertex > &vertices)
Returns face as polygon geometry.
A rectangle specified with double values.
Definition: qgsrectangle.h:41
double y
Definition: qgspoint.h:42
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:486
QPointF toQPointF() const
Returns the point as a QPointF.
Definition: qgspoint.h:264
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
int faceIndexForPoint(const QgsPointXY &point) const
Finds index of triangle at given point It uses spatial indexing.
A class to represent a 2D point.
Definition: qgspointxy.h:43
QList< int > intersects(const QgsRectangle &rectangle) const
Returns a list of face ids with a bounding box which intersects the specified rectangle.
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:111
bool isValid() const
Returns true if the coordinate transform is valid, ie both the source and destination CRS have been s...
QString what() const
Definition: qgsexception.h:48
QVector< QgsMeshVertex > vertices
vertices
const QVector< QgsMeshFace > & triangles() const
Returns triangles.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from...
void setM(double m)
Sets the point&#39;s m-value.
Definition: qgspoint.h:252
CORE_EXPORT std::unique_ptr< QgsPolygon > toPolygon(const QgsMeshFace &face, const QVector< QgsMeshVertex > &vertices)
Returns face as polygon geometry, caller is responsible for delete.
QgsCoordinateTransform coordinateTransform() const
Returns the current coordinate transform for the context.
~QgsTriangularMesh()
Dtor.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
CORE_EXPORT QList< int > nativeFacesFromTriangles(const QList< int > &triangleIndexes, const QVector< int > &trianglesToNativeFaces)
Returns unique native faces indexes from list of triangle indexes.
QList< int > faceIndexesForRectangle(const QgsRectangle &rectangle) const
Finds indexes of triangles intersecting given bounding box It uses spatial indexing.
A spatial index for QgsMeshFace objects.
void update(QgsMesh *nativeMesh, QgsRenderContext *context)
Constructs triangular mesh from layer&#39;s native mesh and context.
Contains information about the context of a rendering operation.
Mesh - vertices and faces.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
QVector< int > QgsMeshFace
List of vertex indexes.
QgsTriangularMesh()
Ctor.
double z
Definition: qgspoint.h:43
QVector< QgsMeshFace > faces
faces
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:65
bool contains(const QgsPointXY *p) const
Tests for containment of a point (uses GEOS)
const QVector< int > & trianglesToNativeFaces() const
Returns mapping between triangles and original faces.
const QVector< QgsMeshVertex > & vertices() const
Returns vertices in map coordinate system.
const QVector< QgsMeshVertex > & centroids() const
Returns centroids of the native faces in map CRS.
double m
Definition: qgspoint.h:44
QgsPoint QgsMeshVertex
xyz coords of vertex
double x
Definition: qgspoint.h:41