QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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 
37 QgsMeshLayerInterpolator::QgsMeshLayerInterpolator( const QgsTriangularMesh &m,
38  const QVector<double> &datasetValues, const QgsMeshDataBlock &activeFaceFlagValues,
39  bool dataIsOnVertices,
40  const QgsRenderContext &context,
41  const QSize &size )
42  : mTriangularMesh( m ),
43  mDatasetValues( datasetValues ),
44  mActiveFaceFlagValues( activeFaceFlagValues ),
45  mContext( context ),
46  mDataOnVertices( dataIsOnVertices ),
47  mOutputSize( size )
48 {
49 }
50 
51 QgsMeshLayerInterpolator::~QgsMeshLayerInterpolator() = default;
52 
53 QgsRasterInterface *QgsMeshLayerInterpolator::clone() const
54 {
55  assert( false ); // we should not need this (hopefully)
56  return nullptr;
57 }
58 
59 Qgis::DataType QgsMeshLayerInterpolator::dataType( int ) const
60 {
61  return Qgis::Float64;
62 }
63 
64 int QgsMeshLayerInterpolator::bandCount() const
65 {
66  return 1;
67 }
68 
69 QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
70 {
71  std::unique_ptr<QgsRasterBlock> outputBlock( new QgsRasterBlock( Qgis::Float64, width, height ) );
72  const double noDataValue = std::numeric_limits<double>::quiet_NaN();
73  outputBlock->setNoDataValue( noDataValue );
74  outputBlock->setIsNoData(); // assume initially that all values are unset
75  double *data = reinterpret_cast<double *>( outputBlock->bits() );
76 
77  const QVector<QgsMeshFace> &triangles = mTriangularMesh.triangles();
78  const QVector<QgsMeshVertex> &vertices = mTriangularMesh.vertices();
79 
80  // currently expecting that triangulation does not add any new extra vertices on the way
81  if ( mDataOnVertices )
82  Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
83 
84  for ( int i = 0; i < triangles.size(); ++i )
85  {
86  if ( feedback && feedback->isCanceled() )
87  break;
88 
89  if ( mContext.renderingStopped() )
90  break;
91 
92  const QgsMeshFace &face = triangles[i];
93 
94  const int v1 = face[0], v2 = face[1], v3 = face[2];
95  const QgsPoint p1 = vertices[v1], p2 = vertices[v2], p3 = vertices[v3];
96 
97  const int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces()[i];
98  const bool isActive = mActiveFaceFlagValues.active( nativeFaceIndex );
99  if ( !isActive )
100  continue;
101 
102  QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox( p1, p2, p3 );
103  if ( !extent.intersects( bbox ) )
104  continue;
105 
106  // Get the BBox of the element in pixels
107  int topLim, bottomLim, leftLim, rightLim;
108  QgsMeshLayerUtils::boundingBoxToScreenRectangle( mContext.mapToPixel(), mOutputSize, bbox, leftLim, rightLim, topLim, bottomLim );
109 
110  // interpolate in the bounding box of the face
111  for ( int j = topLim; j <= bottomLim; j++ )
112  {
113  double *line = data + ( j * width );
114  for ( int k = leftLim; k <= rightLim; k++ )
115  {
116  double val;
117  const QgsPointXY p = mContext.mapToPixel().toMapCoordinates( k, j );
118  if ( mDataOnVertices )
119  val = QgsMeshLayerUtils::interpolateFromVerticesData(
120  p1,
121  p2,
122  p3,
123  mDatasetValues[v1],
124  mDatasetValues[v2],
125  mDatasetValues[v3],
126  p );
127  else
128  {
129  int face = mTriangularMesh.trianglesToNativeFaces()[i];
130  val = QgsMeshLayerUtils::interpolateFromFacesData(
131  p1,
132  p2,
133  p3,
134  mDatasetValues[face],
135  p
136  );
137  }
138 
139  if ( !std::isnan( val ) )
140  {
141  line[k] = val;
142  outputBlock->setIsData( j, k );
143  }
144  }
145  }
146 
147  }
148 
149  return outputBlock.release();
150 }
151 
153 
155  const QgsMeshLayer &layer,
156  const QgsMeshDatasetIndex &datasetIndex,
157  const QgsCoordinateReferenceSystem &destinationCrs,
159  double mapUnitsPerPixel,
160  const QgsRectangle &extent,
161  QgsRasterBlockFeedback *feedback )
162 {
163  if ( !layer.dataProvider() )
164  return nullptr;
165 
166  if ( !datasetIndex.isValid() )
167  return nullptr;
168 
169  int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
170  int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );
171 
172  const QgsPointXY center = extent.center();
173  QgsMapToPixel mapToPixel( mapUnitsPerPixel,
174  center.x(),
175  center.y(),
176  widthPixel,
177  heightPixel,
178  0 );
179  QgsCoordinateTransform transform( layer.crs(), destinationCrs, transformContext );
180 
181  QgsRenderContext renderContext;
182  renderContext.setCoordinateTransform( transform );
183  renderContext.setMapToPixel( mapToPixel );
184  renderContext.setExtent( extent );
185 
186  std::unique_ptr<QgsMesh> nativeMesh = qgis::make_unique<QgsMesh>();
187  layer.dataProvider()->populateMesh( nativeMesh.get() );
188  std::unique_ptr<QgsTriangularMesh> triangularMesh = qgis::make_unique<QgsTriangularMesh>();
189  triangularMesh->update( nativeMesh.get(), &renderContext );
190 
192  bool scalarDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
193 
195  datasetIndex,
196  0,
197  scalarDataOnVertices ? nativeMesh->vertices.count() : nativeMesh->faces.count() );
198 
199  QVector<double> datasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
200  QgsMeshDataBlock activeFaceFlagValues = layer.dataProvider()->areFacesActive(
201  datasetIndex,
202  0,
203  nativeMesh->faces.count() );
204 
205  QgsMeshLayerInterpolator iterpolator(
206  *( triangularMesh.get() ),
207  datasetValues,
208  activeFaceFlagValues,
209  scalarDataOnVertices,
210  renderContext,
211  QSize( widthPixel, heightPixel )
212  );
213 
214  return iterpolator.block( 0, extent, widthPixel, heightPixel, feedback );
215 }
A rectangle specified with double values.
Definition: qgsrectangle.h:41
Triangular/Derived Mesh is mesh with vertices in map coordinates.
QgsCoordinateTransformContext transformContext
Definition: qgsmeshlayer.h:109
void setCoordinateTransform(const QgsCoordinateTransform &t)
Sets the current coordinate transform for the context.
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
QgsMeshDataBlock is a block of integers/doubles that can be used to retrieve: active flags (e...
DataType
Raster data types.
Definition: qgis.h:79
virtual const QgsLayerMetadata & metadata() const
Returns a reference to the layer&#39;s metadata store.
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&#39;s dataset values as raster block.
QgsMesh * nativeMesh()
Returns native mesh (nullptr before rendering)
Sixty four bit floating point (double)
Definition: qgis.h:88
DataType dataType() const
Returns whether dataset group data is defined on vertices or faces.
Perform transforms between map coordinates and device coordinates.
Definition: qgsmaptopixel.h:37
QgsRectangle extent() const override
Returns the extent of the layer.
Raster data container.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
virtual QgsMeshDataBlock datasetValues(QgsMeshDatasetIndex index, int valueIndex, int count) const =0
Returns N vector/scalar values from the index from the dataset.
QgsMeshDataProvider * dataProvider() override
Returns the layer&#39;s data provider, it may be nullptr.
Contains information about the context in which a coordinate transform is executed.
Base class for processing filters like renderers, reprojector, resampler etc.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
virtual void populateMesh(QgsMesh *mesh) const =0
Populates the mesh vertices and faces.
double x
Definition: qgspointxy.h:47
void update(QgsMesh *nativeMesh, QgsRenderContext *context)
Constructs triangular mesh from layer&#39;s native mesh and context.
Contains information about the context of a rendering operation.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
QgsMeshDatasetGroupMetadata is a collection of dataset group metadata such as whether the data is vec...
QVector< int > QgsMeshFace
List of vertex indexes.
This class represents a coordinate reference system (CRS).
QgsMeshDatasetIndex is index that identifies the dataset group (e.g.
Class for doing transforms between two map coordinate systems.
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:90
QgsPointXY center() const
Returns the center point of the rectangle.
Definition: qgsrectangle.h:230
bool intersects(const QgsRectangle &rect) const
Returns true when rectangle intersects with other rectangle.
Definition: qgsrectangle.h:328
Feedback object tailored for raster block reading.
QgsTriangularMesh * triangularMesh()
Returns triangular mesh (nullptr before rendering)
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:85
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
virtual QgsMeshDataBlock areFacesActive(QgsMeshDatasetIndex index, int faceIndex, int count) const =0
Returns whether the faces are active for particular dataset.
virtual QgsMeshDatasetGroupMetadata datasetGroupMetadata(int groupIndex) const =0
Returns dataset group metadata.