QGIS API Documentation  3.23.0-Master (dd0cd13a00)
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 #include "qgsrendercontext.h"
37 
38 QgsMeshLayerInterpolator::QgsMeshLayerInterpolator(
39  const QgsTriangularMesh &m,
40  const QVector<double> &datasetValues,
41  const QgsMeshDataBlock &activeFaceFlagValues,
43  const QgsRenderContext &context,
44  const QSize &size )
45  : mTriangularMesh( m ),
46  mDatasetValues( datasetValues ),
47  mActiveFaceFlagValues( activeFaceFlagValues ),
48  mContext( context ),
49  mDataType( dataType ),
50  mOutputSize( size )
51 {
52 }
53 
54 QgsMeshLayerInterpolator::~QgsMeshLayerInterpolator() = default;
55 
56 QgsRasterInterface *QgsMeshLayerInterpolator::clone() const
57 {
58  assert( false ); // we should not need this (hopefully)
59  return nullptr;
60 }
61 
62 Qgis::DataType QgsMeshLayerInterpolator::dataType( int ) const
63 {
65 }
66 
67 int QgsMeshLayerInterpolator::bandCount() const
68 {
69  return 1;
70 }
71 
72 QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
73 {
74  std::unique_ptr<QgsRasterBlock> outputBlock( new QgsRasterBlock( Qgis::DataType::Float64, width, height ) );
75  const double noDataValue = std::numeric_limits<double>::quiet_NaN();
76  outputBlock->setNoDataValue( noDataValue );
77  outputBlock->setIsNoData(); // assume initially that all values are unset
78  double *data = reinterpret_cast<double *>( outputBlock->bits() );
79 
80  QList<int> spatialIndexTriangles;
81  int indexCount;
82  if ( mSpatialIndexActive )
83  {
84  spatialIndexTriangles = mTriangularMesh.faceIndexesForRectangle( extent );
85  indexCount = spatialIndexTriangles.count();
86  }
87  else
88  {
89  indexCount = mTriangularMesh.triangles().count();
90  }
91 
92  if ( mTriangularMesh.contains( QgsMesh::ElementType::Edge ) )
93  {
94  return outputBlock.release();
95  }
96 
97  const QVector<QgsMeshVertex> &vertices = mTriangularMesh.vertices();
98 
99  // currently expecting that triangulation does not add any new extra vertices on the way
100  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
101  Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
102 
103  for ( int i = 0; i < indexCount; ++i )
104  {
105  if ( feedback && feedback->isCanceled() )
106  break;
107 
108  if ( mContext.renderingStopped() )
109  break;
110 
111  int triangleIndex;
112  if ( mSpatialIndexActive )
113  triangleIndex = spatialIndexTriangles[i];
114  else
115  triangleIndex = i;
116 
117  const QgsMeshFace &face = mTriangularMesh.triangles()[triangleIndex];
118 
119  if ( face.isEmpty() )
120  continue;
121 
122  const int v1 = face[0], v2 = face[1], v3 = face[2];
123  const QgsPointXY &p1 = vertices[v1], &p2 = vertices[v2], &p3 = vertices[v3];
124 
125  const int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
126  const bool isActive = mActiveFaceFlagValues.active( nativeFaceIndex );
127  if ( !isActive )
128  continue;
129 
130  const QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox( p1, p2, p3 );
131  if ( !extent.intersects( bbox ) )
132  continue;
133 
134  // Get the BBox of the element in pixels
135  int topLim, bottomLim, leftLim, rightLim;
136  QgsMeshLayerUtils::boundingBoxToScreenRectangle( mContext.mapToPixel(), mOutputSize, bbox, leftLim, rightLim, topLim, bottomLim );
137 
138  double value( 0 ), value1( 0 ), value2( 0 ), value3( 0 );
139  const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
140 
141  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
142  {
143  value1 = mDatasetValues[v1];
144  value2 = mDatasetValues[v2];
145  value3 = mDatasetValues[v3];
146  }
147  else
148  value = mDatasetValues[faceIdx];
149 
150  // interpolate in the bounding box of the face
151  for ( int j = topLim; j <= bottomLim; j++ )
152  {
153  double *line = data + ( j * width );
154  for ( int k = leftLim; k <= rightLim; k++ )
155  {
156  double val;
157  const QgsPointXY p = mContext.mapToPixel().toMapCoordinates( k, j );
158  if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
159  val = QgsMeshLayerUtils::interpolateFromVerticesData(
160  p1,
161  p2,
162  p3,
163  value1,
164  value2,
165  value3,
166  p );
167  else
168  {
169  val = QgsMeshLayerUtils::interpolateFromFacesData(
170  p1,
171  p2,
172  p3,
173  value,
174  p
175  );
176  }
177  if ( !std::isnan( val ) )
178  {
179  line[k] = val;
180  outputBlock->setIsData( j, k );
181  }
182  }
183  }
184 
185  }
186 
187  return outputBlock.release();
188 }
189 
190 void QgsMeshLayerInterpolator::setSpatialIndexActive( bool active ) {mSpatialIndexActive = active;}
191 
193 
195  const QgsMeshLayer &layer,
196  const QgsMeshDatasetIndex &datasetIndex,
197  const QgsCoordinateReferenceSystem &destinationCrs,
198  const QgsCoordinateTransformContext &transformContext,
199  double mapUnitsPerPixel,
200  const QgsRectangle &extent,
201  QgsRasterBlockFeedback *feedback )
202 {
203  if ( !layer.dataProvider() )
204  return nullptr;
205 
206  if ( !datasetIndex.isValid() )
207  return nullptr;
208 
209  const int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
210  const int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
211 
212  const QgsPointXY center = extent.center();
213  const QgsMapToPixel mapToPixel( mapUnitsPerPixel,
214  center.x(),
215  center.y(),
216  widthPixel,
217  heightPixel,
218  0 );
219  const QgsCoordinateTransform transform( layer.crs(), destinationCrs, transformContext );
220 
221  QgsRenderContext renderContext;
222  renderContext.setCoordinateTransform( transform );
223  renderContext.setMapToPixel( mapToPixel );
224  renderContext.setExtent( extent );
225 
226  std::unique_ptr<QgsMesh> nativeMesh = std::make_unique<QgsMesh>();
227  layer.dataProvider()->populateMesh( nativeMesh.get() );
228  std::unique_ptr<QgsTriangularMesh> triangularMesh = std::make_unique<QgsTriangularMesh>();
229  triangularMesh->update( nativeMesh.get(), transform );
230 
231  const QgsMeshDatasetGroupMetadata metadata = layer.dataProvider()->datasetGroupMetadata( datasetIndex );
232  const QgsMeshDatasetGroupMetadata::DataType scalarDataType = QgsMeshLayerUtils::datasetValuesType( metadata.dataType() );
233  const int count = QgsMeshLayerUtils::datasetValuesCount( nativeMesh.get(), scalarDataType );
234  const QgsMeshDataBlock vals = QgsMeshLayerUtils::datasetValues(
235  &layer,
236  datasetIndex,
237  0,
238  count );
239  if ( !vals.isValid() )
240  return nullptr;
241 
242  const QVector<double> datasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
243  const QgsMeshDataBlock activeFaceFlagValues = layer.dataProvider()->areFacesActive(
244  datasetIndex,
245  0,
246  nativeMesh->faces.count() );
247 
248  QgsMeshLayerInterpolator interpolator(
249  *( triangularMesh.get() ),
250  datasetValues,
251  activeFaceFlagValues,
252  scalarDataType,
253  renderContext,
254  QSize( widthPixel, heightPixel )
255  );
256 
257  return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
258 }
259 
261  const QgsTriangularMesh &triangularMesh,
262  const QgsMeshDataBlock &datasetValues,
263  const QgsMeshDataBlock &activeFlags,
265  const QgsCoordinateTransform &transform,
266  double mapUnitsPerPixel,
267  const QgsRectangle &extent,
268  QgsRasterBlockFeedback *feedback )
269 {
270 
271  const int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
272  const int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
273 
274  const QgsPointXY center = extent.center();
275  const QgsMapToPixel mapToPixel( mapUnitsPerPixel,
276  center.x(),
277  center.y(),
278  widthPixel,
279  heightPixel,
280  0 );
281 
282  QgsRenderContext renderContext;
283  renderContext.setCoordinateTransform( transform );
284  renderContext.setMapToPixel( mapToPixel );
285  renderContext.setExtent( extent );
286 
287  const QVector<double> magnitudes = QgsMeshLayerUtils::calculateMagnitudes( datasetValues );
288 
289  QgsMeshLayerInterpolator interpolator(
290  triangularMesh,
291  magnitudes,
292  activeFlags,
293  dataType,
294  renderContext,
295  QSize( widthPixel, heightPixel )
296  );
297 
298  return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
299 }
DataType
Raster data types.
Definition: qgis.h:121
@ Float64
Sixty four bit floating point (double)
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 SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
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:97
QgsMeshDataProvider * dataProvider() override
Returns the layer's data provider, it may be nullptr.
A class to represent a 2D point.
Definition: qgspointxy.h:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
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 SIP_HOLDGIL
Returns true when rectangle intersects with other rectangle.
Definition: qgsrectangle.h:349
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
QgsPointXY center() const SIP_HOLDGIL
Returns the center point of the rectangle.
Definition: qgsrectangle.h:251
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.