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