QGIS API Documentation 3.43.0-Master (b60ef06885e)
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 "qgspointxy.h"
29#include "qgsmeshlayerutils.h"
30#include "qgsmeshlayer.h"
32#include "qgsmeshdataprovider.h"
33#include "qgsrendercontext.h"
34#include "qgselevationmap.h"
35
36QgsMeshLayerInterpolator::QgsMeshLayerInterpolator(
37 const QgsTriangularMesh &m,
38 const QVector<double> &datasetValues,
39 const QgsMeshDataBlock &activeFaceFlagValues,
41 const QgsRenderContext &context,
42 const QSize &size )
43 : mTriangularMesh( m ),
44 mDatasetValues( datasetValues ),
45 mActiveFaceFlagValues( activeFaceFlagValues ),
46 mContext( context ),
47 mDataType( dataType ),
48 mOutputSize( size )
49{
50}
51
52QgsMeshLayerInterpolator::~QgsMeshLayerInterpolator() = default;
53
54QgsRasterInterface *QgsMeshLayerInterpolator::clone() const
55{
56 assert( false ); // we should not need this (hopefully)
57 return nullptr;
58}
59
60Qgis::DataType QgsMeshLayerInterpolator::dataType( int ) const
61{
63}
64
65int QgsMeshLayerInterpolator::bandCount() const
66{
67 return 1;
68}
69
70QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
71{
72 auto outputBlock = std::make_unique<QgsRasterBlock>( Qgis::DataType::Float64, width, height );
73 const double noDataValue = std::numeric_limits<double>::quiet_NaN();
74 outputBlock->setNoDataValue( noDataValue );
75 outputBlock->setIsNoData(); // assume initially that all values are unset
76 double *data = reinterpret_cast<double *>( outputBlock->bits() );
77
78 QList<int> spatialIndexTriangles;
79 int indexCount;
80 if ( mSpatialIndexActive )
81 {
82 spatialIndexTriangles = mTriangularMesh.faceIndexesForRectangle( extent );
83 indexCount = spatialIndexTriangles.count();
84 }
85 else
86 {
87 indexCount = mTriangularMesh.triangles().count();
88 }
89
90 if ( mTriangularMesh.contains( QgsMesh::ElementType::Edge ) )
91 {
92 return outputBlock.release();
93 }
94
95 const QVector<QgsMeshVertex> &vertices = mTriangularMesh.vertices();
96
97 // currently expecting that triangulation does not add any new extra vertices on the way
99 Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
100
101 double pixelRatio = mContext.devicePixelRatio();
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, pixelRatio );
137
138 double value( 0 ), value1( 0 ), value2( 0 ), value3( 0 );
139 const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
140
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 / pixelRatio, j / pixelRatio );
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 if ( mRenderElevation )
187 {
188 QgsElevationMap *elevationMap = mContext.elevationMap();
189 if ( elevationMap && elevationMap->isValid() )
190 elevationMap->fillWithRasterBlock( outputBlock.get(), 0, 0, mElevationScale, mElevationOffset );
191 }
192
193 return outputBlock.release();
194}
195
196void QgsMeshLayerInterpolator::setSpatialIndexActive( bool active )
197{
198 mSpatialIndexActive = active;
199}
200
201void QgsMeshLayerInterpolator::setElevationMapSettings( bool renderElevationMap, double elevationScale, double elevationOffset )
202{
203 mRenderElevation = renderElevationMap;
204 mElevationScale = elevationScale;
205 mElevationOffset = elevationOffset;
206}
207
DataType
Raster data types.
Definition qgis.h:351
@ Float64
Sixty four bit floating point (double)
Stores a digital elevation model in a raster image which may get updated as a part of the map layer r...
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
A block of integers/doubles from a mesh dataset.
DataType
Location of where data is specified for datasets in the dataset group.
@ DataOnVertices
Data is defined on vertices.
Represents a 2D point.
Definition qgspointxy.h:60
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.
Contains information about the context of a rendering operation.
A triangular/derived mesh with vertices in map coordinates.
QVector< int > QgsMeshFace
List of vertex indexes.