QGIS API Documentation  3.18.1-Zürich (202f1bf7e5)
qgsmeshlayerinterpolator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsmeshlayerinterpolator.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 
19 
20 #include <memory>
21 #include <limits>
22 
24 
25 #include "qgis.h"
26 #include "qgsrasterinterface.h"
27 #include "qgsmaptopixel.h"
28 #include "qgsvector.h"
29 #include "qgspoint.h"
30 #include "qgspointxy.h"
31 #include "qgsmeshlayerutils.h"
32 #include "qgsmeshlayer.h"
34 #include "qgscoordinatetransform.h"
35 #include "qgsmeshdataprovider.h"
36 
37 QgsMeshLayerInterpolator::QgsMeshLayerInterpolator(
38  const QgsTriangularMesh &m,
39  const QVector<double> &datasetValues,
40  const QgsMeshDataBlock &activeFaceFlagValues,
42  const QgsRenderContext &context,
43  const QSize &size )
44  : mTriangularMesh( m ),
45  mDatasetValues( datasetValues ),
46  mActiveFaceFlagValues( activeFaceFlagValues ),
47  mContext( context ),
48  mDataType( dataType ),
49  mOutputSize( size )
50 {
51 }
52 
53 QgsMeshLayerInterpolator::~QgsMeshLayerInterpolator() = default;
54 
55 QgsRasterInterface *QgsMeshLayerInterpolator::clone() const
56 {
57  assert( false ); // we should not need this (hopefully)
58  return nullptr;
59 }
60 
61 Qgis::DataType QgsMeshLayerInterpolator::dataType( int ) const
62 {
63  return Qgis::Float64;
64 }
65 
66 int QgsMeshLayerInterpolator::bandCount() const
67 {
68  return 1;
69 }
70 
71 QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
72 {
73  std::unique_ptr<QgsRasterBlock> outputBlock( new QgsRasterBlock( Qgis::Float64, width, height ) );
74  const double noDataValue = std::numeric_limits<double>::quiet_NaN();
75  outputBlock->setNoDataValue( noDataValue );
76  outputBlock->setIsNoData(); // assume initially that all values are unset
77  double *data = reinterpret_cast<double *>( outputBlock->bits() );
78 
79  QList<int> spatialIndexTriangles;
80  int indexCount;
81  if ( mSpatialIndexActive )
82  {
83  spatialIndexTriangles = mTriangularMesh.faceIndexesForRectangle( extent );
84  indexCount = spatialIndexTriangles.count();
85  }
86  else
87  {
88  indexCount = mTriangularMesh.triangles().count();
89  }
90 
91  if ( mTriangularMesh.contains( QgsMesh::ElementType::Edge ) )
92  {
93  return outputBlock.release();
94  }
95 
96  const QVector<QgsMeshVertex> &vertices = mTriangularMesh.vertices();
97 
98  // currently expecting that triangulation does not add any new extra vertices on the way
99  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
100  Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
101 
102  for ( int i = 0; i < indexCount; ++i )
103  {
104  if ( feedback && feedback->isCanceled() )
105  break;
106 
107  if ( mContext.renderingStopped() )
108  break;
109 
110  int triangleIndex;
111  if ( mSpatialIndexActive )
112  triangleIndex = spatialIndexTriangles[i];
113  else
114  triangleIndex = i;
115 
116  const QgsMeshFace &face = mTriangularMesh.triangles()[triangleIndex];
117 
118  const int v1 = face[0], v2 = face[1], v3 = face[2];
119  const QgsPointXY &p1 = vertices[v1], &p2 = vertices[v2], &p3 = vertices[v3];
120 
121  const int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
122  const bool isActive = mActiveFaceFlagValues.active( nativeFaceIndex );
123  if ( !isActive )
124  continue;
125 
126  QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox( p1, p2, p3 );
127  if ( !extent.intersects( bbox ) )
128  continue;
129 
130  // Get the BBox of the element in pixels
131  int topLim, bottomLim, leftLim, rightLim;
132  QgsMeshLayerUtils::boundingBoxToScreenRectangle( mContext.mapToPixel(), mOutputSize, bbox, leftLim, rightLim, topLim, bottomLim );
133 
134  double value( 0 ), value1( 0 ), value2( 0 ), value3( 0 );
135  const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
136 
137  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
138  {
139  value1 = mDatasetValues[v1];
140  value2 = mDatasetValues[v2];
141  value3 = mDatasetValues[v3];
142  }
143  else
144  value = mDatasetValues[faceIdx];
145 
146  // interpolate in the bounding box of the face
147  for ( int j = topLim; j <= bottomLim; j++ )
148  {
149  double *line = data + ( j * width );
150  for ( int k = leftLim; k <= rightLim; k++ )
151  {
152  double val;
153  const QgsPointXY p = mContext.mapToPixel().toMapCoordinates( k, j );
154  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
155  val = QgsMeshLayerUtils::interpolateFromVerticesData(
156  p1,
157  p2,
158  p3,
159  value1,
160  value2,
161  value3,
162  p );
163  else
164  {
165  val = QgsMeshLayerUtils::interpolateFromFacesData(
166  p1,
167  p2,
168  p3,
169  value,
170  p
171  );
172  }
173  if ( !std::isnan( val ) )
174  {
175  line[k] = val;
176  outputBlock->setIsData( j, k );
177  }
178  }
179  }
180 
181  }
182 
183  return outputBlock.release();
184 }
185 
186 void QgsMeshLayerInterpolator::setSpatialIndexActive( bool active ) {mSpatialIndexActive = active;}
187 
189 
191  const QgsMeshLayer &layer,
192  const QgsMeshDatasetIndex &datasetIndex,
193  const QgsCoordinateReferenceSystem &destinationCrs,
194  const QgsCoordinateTransformContext &transformContext,
195  double mapUnitsPerPixel,
196  const QgsRectangle &extent,
197  QgsRasterBlockFeedback *feedback )
198 {
199  if ( !layer.dataProvider() )
200  return nullptr;
201 
202  if ( !datasetIndex.isValid() )
203  return nullptr;
204 
205  int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
206  int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
207 
208  const QgsPointXY center = extent.center();
209  QgsMapToPixel mapToPixel( mapUnitsPerPixel,
210  center.x(),
211  center.y(),
212  widthPixel,
213  heightPixel,
214  0 );
215  QgsCoordinateTransform transform( layer.crs(), destinationCrs, transformContext );
216 
217  QgsRenderContext renderContext;
218  renderContext.setCoordinateTransform( transform );
219  renderContext.setMapToPixel( mapToPixel );
220  renderContext.setExtent( extent );
221 
222  std::unique_ptr<QgsMesh> nativeMesh = qgis::make_unique<QgsMesh>();
223  layer.dataProvider()->populateMesh( nativeMesh.get() );
224  std::unique_ptr<QgsTriangularMesh> triangularMesh = qgis::make_unique<QgsTriangularMesh>();
225  triangularMesh->update( nativeMesh.get(), transform );
226 
227  const QgsMeshDatasetGroupMetadata metadata = layer.dataProvider()->datasetGroupMetadata( datasetIndex );
228  QgsMeshDatasetGroupMetadata::DataType scalarDataType = QgsMeshLayerUtils::datasetValuesType( metadata.dataType() );
229  const int count = QgsMeshLayerUtils::datasetValuesCount( nativeMesh.get(), scalarDataType );
230  QgsMeshDataBlock vals = QgsMeshLayerUtils::datasetValues(
231  &layer,
232  datasetIndex,
233  0,
234  count );
235  if ( !vals.isValid() )
236  return nullptr;
237 
238  QVector<double> datasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
239  QgsMeshDataBlock activeFaceFlagValues = layer.dataProvider()->areFacesActive(
240  datasetIndex,
241  0,
242  nativeMesh->faces.count() );
243 
244  QgsMeshLayerInterpolator interpolator(
245  *( triangularMesh.get() ),
246  datasetValues,
247  activeFaceFlagValues,
248  scalarDataType,
249  renderContext,
250  QSize( widthPixel, heightPixel )
251  );
252 
253  return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
254 }
255 
257  const QgsTriangularMesh &triangularMesh,
258  const QgsMeshDataBlock &datasetValues,
259  const QgsMeshDataBlock &activeFlags,
261  const QgsCoordinateTransform &transform,
262  double mapUnitsPerPixel,
263  const QgsRectangle &extent,
264  QgsRasterBlockFeedback *feedback )
265 {
266 
267  int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
268  int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
269 
270  const QgsPointXY center = extent.center();
271  QgsMapToPixel mapToPixel( mapUnitsPerPixel,
272  center.x(),
273  center.y(),
274  widthPixel,
275  heightPixel,
276  0 );
277 
278  QgsRenderContext renderContext;
279  renderContext.setCoordinateTransform( transform );
280  renderContext.setMapToPixel( mapToPixel );
281  renderContext.setExtent( extent );
282 
283  QVector<double> magnitudes = QgsMeshLayerUtils::calculateMagnitudes( datasetValues );
284 
285  QgsMeshLayerInterpolator interpolator(
286  triangularMesh,
287  magnitudes,
288  activeFlags,
289  dataType,
290  renderContext,
291  QSize( widthPixel, heightPixel )
292  );
293 
294  return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
295 }
DataType
Raster data types.
Definition: qgis.h:102
@ Float64
Sixty four bit floating point (double)
Definition: qgis.h:110
This class represents a coordinate reference system (CRS).
Contains information about the context in which a coordinate transform is executed.
Class for doing transforms between two map coordinate systems.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:91
Perform transforms between map coordinates and device coordinates.
Definition: qgsmaptopixel.h:39
QgsMeshDataBlock is a block of integers/doubles that can be used to retrieve: active flags (e....
bool isValid() const
Whether the block is valid.
virtual void populateMesh(QgsMesh *mesh) const =0
Populates the mesh vertices, edges and faces.
QgsMeshDatasetGroupMetadata is a collection of dataset group metadata such as whether the data is vec...
DataType dataType() const
Returns whether dataset group data is defined on vertices or faces or volumes.
DataType
Location of where data is specified for datasets in the dataset group.
QgsMeshDatasetIndex is index that identifies the dataset group (e.g.
bool isValid() const
Returns whether index is valid, ie at least groups is set.
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:95
QgsMeshDataProvider * dataProvider() override
Returns the layer's data provider, it may be nullptr.
A class to represent a 2D point.
Definition: qgspointxy.h:44
double y
Definition: qgspointxy.h:48
Q_GADGET double x
Definition: qgspointxy.h:47
Feedback object tailored for raster block reading.
Raster data container.
Base class for processing filters like renderers, reprojector, resampler etc.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
bool intersects(const QgsRectangle &rect) const
Returns true when rectangle intersects with other rectangle.
Definition: qgsrectangle.h:328
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
QgsPointXY center() const SIP_HOLDGIL
Returns the center point of the rectangle.
Definition: qgsrectangle.h:230
Contains information about the context of a rendering operation.
void setCoordinateTransform(const QgsCoordinateTransform &t)
Sets the current coordinate transform for the context.
void setExtent(const QgsRectangle &extent)
When rendering a map layer, calling this method sets the "clipping" extent for the layer (in the laye...
void setMapToPixel(const QgsMapToPixel &mtp)
Sets the context's map to pixel transform, which transforms between map coordinates and device coordi...
Triangular/Derived Mesh is mesh with vertices in map coordinates.
CORE_EXPORT QgsRasterBlock * exportRasterBlock(const QgsMeshLayer &layer, const QgsMeshDatasetIndex &datasetIndex, const QgsCoordinateReferenceSystem &destinationCrs, const QgsCoordinateTransformContext &transformContext, double mapUnitsPerPixel, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Exports mesh layer's dataset values as raster block.
QVector< int > QgsMeshFace
List of vertex indexes.