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