QGIS API Documentation 3.37.0-Master (c2d15952569)
Loading...
Searching...
No Matches
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
102 Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
103
104 double pixelRatio = mContext.devicePixelRatio();
105
106 for ( int i = 0; i < indexCount; ++i )
107 {
108 if ( feedback && feedback->isCanceled() )
109 break;
110
111 if ( mContext.renderingStopped() )
112 break;
113
114 int triangleIndex;
115 if ( mSpatialIndexActive )
116 triangleIndex = spatialIndexTriangles[i];
117 else
118 triangleIndex = i;
119
120 const QgsMeshFace &face = mTriangularMesh.triangles()[triangleIndex];
121
122 if ( face.isEmpty() )
123 continue;
124
125 const int v1 = face[0], v2 = face[1], v3 = face[2];
126 const QgsPointXY &p1 = vertices[v1], &p2 = vertices[v2], &p3 = vertices[v3];
127
128 const int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
129 const bool isActive = mActiveFaceFlagValues.active( nativeFaceIndex );
130 if ( !isActive )
131 continue;
132
133 const QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox( p1, p2, p3 );
134 if ( !extent.intersects( bbox ) )
135 continue;
136
137 // Get the BBox of the element in pixels
138 int topLim, bottomLim, leftLim, rightLim;
139 QgsMeshLayerUtils::boundingBoxToScreenRectangle( mContext.mapToPixel(), mOutputSize, bbox, leftLim, rightLim, topLim, bottomLim, pixelRatio );
140
141 double value( 0 ), value1( 0 ), value2( 0 ), value3( 0 );
142 const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
143
145 {
146 value1 = mDatasetValues[v1];
147 value2 = mDatasetValues[v2];
148 value3 = mDatasetValues[v3];
149 }
150 else
151 value = mDatasetValues[faceIdx];
152
153 // interpolate in the bounding box of the face
154 for ( int j = topLim; j <= bottomLim; j++ )
155 {
156 double *line = data + ( j * width );
157 for ( int k = leftLim; k <= rightLim; k++ )
158 {
159 double val;
160 const QgsPointXY p = mContext.mapToPixel().toMapCoordinates( k / pixelRatio, j / pixelRatio );
162 val = QgsMeshLayerUtils::interpolateFromVerticesData(
163 p1,
164 p2,
165 p3,
166 value1,
167 value2,
168 value3,
169 p );
170 else
171 {
172 val = QgsMeshLayerUtils::interpolateFromFacesData(
173 p1,
174 p2,
175 p3,
176 value,
177 p
178 );
179 }
180 if ( !std::isnan( val ) )
181 {
182 line[k] = val;
183 outputBlock->setIsData( j, k );
184 }
185 }
186 }
187 }
188
189 if ( mRenderElevation )
190 {
191 QgsElevationMap *elevationMap = mContext.elevationMap();
192 if ( elevationMap && elevationMap->isValid() )
193 elevationMap->fillWithRasterBlock( outputBlock.get(), 0, 0, mElevationScale, mElevationOffset );
194 }
195
196 return outputBlock.release();
197}
198
199void QgsMeshLayerInterpolator::setSpatialIndexActive( bool active )
200{
201 mSpatialIndexActive = active;
202}
203
204void QgsMeshLayerInterpolator::setElevationMapSettings( bool renderElevationMap, double elevationScale, double elevationOffset )
205{
206 mRenderElevation = renderElevationMap;
207 mElevationScale = elevationScale;
208 mElevationOffset = elevationOffset;
209}
210
212
214 const QgsMeshLayer &layer,
215 const QgsMeshDatasetIndex &datasetIndex,
216 const QgsCoordinateReferenceSystem &destinationCrs,
217 const QgsCoordinateTransformContext &transformContext,
218 double mapUnitsPerPixel,
219 const QgsRectangle &extent,
220 QgsRasterBlockFeedback *feedback )
221{
222 if ( !layer.dataProvider() )
223 return nullptr;
224
225 if ( !datasetIndex.isValid() )
226 return nullptr;
227
228 const int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
229 const int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
230
231 const QgsPointXY center = extent.center();
232 const QgsMapToPixel mapToPixel( mapUnitsPerPixel,
233 center.x(),
234 center.y(),
235 widthPixel,
236 heightPixel,
237 0 );
238 const QgsCoordinateTransform transform( layer.crs(), destinationCrs, transformContext );
239
240 QgsRenderContext renderContext;
241 renderContext.setCoordinateTransform( transform );
242 renderContext.setMapToPixel( mapToPixel );
243 renderContext.setExtent( extent );
244
245 std::unique_ptr<QgsMesh> nativeMesh = std::make_unique<QgsMesh>();
246 layer.dataProvider()->populateMesh( nativeMesh.get() );
247 std::unique_ptr<QgsTriangularMesh> triangularMesh = std::make_unique<QgsTriangularMesh>();
248 triangularMesh->update( nativeMesh.get(), transform );
249
250 const QgsMeshDatasetGroupMetadata metadata = layer.datasetGroupMetadata( datasetIndex );
251 const QgsMeshDatasetGroupMetadata::DataType scalarDataType = QgsMeshLayerUtils::datasetValuesType( metadata.dataType() );
252 const int count = QgsMeshLayerUtils::datasetValuesCount( nativeMesh.get(), scalarDataType );
253 const QgsMeshDataBlock vals = QgsMeshLayerUtils::datasetValues(
254 &layer,
255 datasetIndex,
256 0,
257 count );
258 if ( !vals.isValid() )
259 return nullptr;
260
261 const QVector<double> datasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
262 const QgsMeshDataBlock activeFaceFlagValues = layer.areFacesActive(
263 datasetIndex,
264 0,
265 nativeMesh->faces.count() );
266
267 QgsMeshLayerInterpolator interpolator(
268 *( triangularMesh.get() ),
269 datasetValues,
270 activeFaceFlagValues,
271 scalarDataType,
272 renderContext,
273 QSize( widthPixel, heightPixel )
274 );
275
276 return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
277}
278
280 const QgsTriangularMesh &triangularMesh,
281 const QgsMeshDataBlock &datasetValues,
282 const QgsMeshDataBlock &activeFlags,
284 const QgsCoordinateTransform &transform,
285 double mapUnitsPerPixel,
286 const QgsRectangle &extent,
287 QgsRasterBlockFeedback *feedback )
288{
289
290 const int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
291 const int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
292
293 const QgsPointXY center = extent.center();
294 const QgsMapToPixel mapToPixel( mapUnitsPerPixel,
295 center.x(),
296 center.y(),
297 widthPixel,
298 heightPixel,
299 0 );
300
301 QgsRenderContext renderContext;
302 renderContext.setCoordinateTransform( transform );
303 renderContext.setMapToPixel( mapToPixel );
304 renderContext.setExtent( extent );
305
306 const QVector<double> magnitudes = QgsMeshLayerUtils::calculateMagnitudes( datasetValues );
307
308 QgsMeshLayerInterpolator interpolator(
309 triangularMesh,
310 magnitudes,
311 activeFlags,
312 dataType,
313 renderContext,
314 QSize( widthPixel, heightPixel )
315 );
316
317 return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
318}
DataType
Raster data types.
Definition qgis.h:269
@ 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
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:81
Perform transforms between map coordinates and device coordinates.
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.
@ DataOnVertices
Data is defined on vertices.
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.
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:60
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
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.
bool intersects(const QgsRectangle &rect) const
Returns true when rectangle intersects with other rectangle.
double width() const
Returns the width of the rectangle.
QgsPointXY center() const
Returns the center point of the rectangle.
double height() const
Returns the height of the rectangle.
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.