[mesh] [feature] add function to identify value on the point

This commit is contained in:
Peter Petrik 2018-08-09 10:40:49 +02:00
parent 51b63e6b46
commit c79e1d0601
11 changed files with 174 additions and 18 deletions

View File

@ -181,6 +181,11 @@ Returns whether dataset group has scalar data
bool isOnVertices() const;
%Docstring
Returns whether dataset group data is defined on vertices
%End
bool isOnFaces() const;
%Docstring
Returns whether dataset group data is defined on faces
%End
};

View File

@ -189,6 +189,30 @@ If dataset is not vector based, do nothing. Triggers repaint
QgsMeshDatasetIndex activeVectorDataset() const;
%Docstring
Returns active vector dataset
%End
QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;
%Docstring
Interpolates the value on the given point from given dataset.
Returns NaN for NaN values, for values outside range or when
triangular mesh is not yet initialized on given point
%End
signals:
void activeScalarDatasetChanged( QgsMeshDatasetIndex index );
%Docstring
Emitted when active scalar dataset is changed
.. versionadded:: 3.4
%End
void activeVectorDatasetChanged( QgsMeshDatasetIndex index );
%Docstring
Emitted when active vector dataset is changed
.. versionadded:: 3.4
%End
private: // Private methods

View File

@ -172,6 +172,11 @@ bool QgsMeshDatasetGroupMetadata::isOnVertices() const
return mIsOnVertices;
}
bool QgsMeshDatasetGroupMetadata::isOnFaces() const
{
return !mIsOnVertices;
}
int QgsMeshDatasetSourceInterface::datasetCount( QgsMeshDatasetIndex index ) const
{
return datasetCount( index.group() );

View File

@ -168,6 +168,11 @@ class CORE_EXPORT QgsMeshDatasetGroupMetadata
*/
bool isOnVertices() const;
/**
* \brief Returns whether dataset group data is defined on faces
*/
bool isOnFaces() const;
private:
QString mName;
bool mIsScalar = false;

View File

@ -16,6 +16,7 @@
***************************************************************************/
#include <cstddef>
#include <limits>
#include <QUuid>
@ -27,6 +28,7 @@
#include "qgsproviderregistry.h"
#include "qgsreadwritecontext.h"
#include "qgstriangularmesh.h"
#include "qgsmeshlayerinterpolator.h"
QgsMeshLayer::QgsMeshLayer( const QString &meshLayerPath,
const QString &baseName,
@ -154,6 +156,8 @@ void QgsMeshLayer::setActiveScalarDataset( QgsMeshDatasetIndex index )
mActiveScalarDataset = QgsMeshDatasetIndex();
triggerRepaint();
emit activeScalarDatasetChanged( mActiveScalarDataset );
}
void QgsMeshLayer::setActiveVectorDataset( QgsMeshDatasetIndex index )
@ -175,6 +179,44 @@ void QgsMeshLayer::setActiveVectorDataset( QgsMeshDatasetIndex index )
}
triggerRepaint();
emit activeVectorDatasetChanged( mActiveVectorDataset );
}
QgsMeshDatasetValue QgsMeshLayer::datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const
{
QgsMeshDatasetValue value;
if ( mTriangularMesh && dataProvider() && dataProvider()->isValid() && index.isValid() )
{
int face_index = mTriangularMesh->faceIndexForPoint( point ) ;
if ( face_index >= 0 )
{
bool isOnFaces = dataProvider()->datasetGroupMetadata( index ).isOnFaces();
if ( isOnFaces )
{
int native_face_index = mTriangularMesh->trianglesToNativeFaces().at( face_index );
return dataProvider()->datasetValue( index, native_face_index );
}
else
{
const QgsMeshFace &face = mTriangularMesh->triangles()[face_index];
const int v1 = face[0], v2 = face[1], v3 = face[2];
const QgsPoint p1 = mTriangularMesh->vertices()[v1], p2 = mTriangularMesh->vertices()[v2], p3 = mTriangularMesh->vertices()[v3];
const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
const QgsMeshDatasetValue val3 = dataProvider()->datasetValue( index, v3 );
const double x = QgsMeshLayerInterpolator::interpolateFromVerticesData( p1, p2, p3, val1.x(), val2.x(), val3.x(), point );
double y = std::numeric_limits<double>::quiet_NaN();
bool isVector = dataProvider()->datasetGroupMetadata( index ).isVector();
if ( isVector )
y = QgsMeshLayerInterpolator::interpolateFromVerticesData( p1, p2, p3, val1.y(), val2.y(), val3.y(), point );
return QgsMeshDatasetValue( x, y );
}
}
}
return value;
}
void QgsMeshLayer::fillNativeMesh()

View File

@ -183,6 +183,30 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
//! Returns active vector dataset
QgsMeshDatasetIndex activeVectorDataset() const { return mActiveVectorDataset; }
/**
* Interpolates the value on the given point from given dataset.
*
* Returns NaN for NaN values, for values outside range or when
* triangular mesh is not yet initialized on given point
*/
QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;
signals:
/**
* Emitted when active scalar dataset is changed
*
* \since QGIS 3.4
*/
void activeScalarDatasetChanged( QgsMeshDatasetIndex index );
/**
* Emitted when active vector dataset is changed
*
* \since QGIS 3.4
*/
void activeVectorDatasetChanged( QgsMeshDatasetIndex index );
private: // Private methods
/**

View File

@ -88,9 +88,8 @@ static bool E3T_physicalToBarycentric( const QgsPointXY &pA, const QgsPointXY &p
return true;
}
double interpolateFromVerticesData( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
double val1, double val2, double val3, const QgsPointXY &pt )
double QgsMeshLayerInterpolator::interpolateFromVerticesData( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
double val1, double val2, double val3, const QgsPointXY &pt )
{
double lam1, lam2, lam3;
if ( !E3T_physicalToBarycentric( p1, p2, p3, pt, lam1, lam2, lam3 ) )

View File

@ -59,6 +59,12 @@ class QgsMeshLayerInterpolator : public QgsRasterInterface
int bandCount() const override;
QgsRasterBlock *block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback = nullptr ) override;
static double interpolateFromVerticesData( const QgsPointXY &p1,
const QgsPointXY &p2,
const QgsPointXY &p3,
double val1, double val2, double val3, const QgsPointXY &pt );
private:
const QgsTriangularMesh &mTriangularMesh;
const QVector<double> &mDatasetValues;

View File

@ -229,17 +229,7 @@ void QgsMeshLayerRenderer::renderMesh( const std::unique_ptr<QgsSymbol> &symbol,
const QgsMeshFace &face = faces[i];
QgsFeature feat;
feat.setFields( fields );
QVector<QgsPointXY> ring;
for ( int j = 0; j < face.size(); ++j )
{
int vertex_id = face[j];
Q_ASSERT( vertex_id < mTriangularMesh.vertices().size() ); //Triangular mesh vertices contains also native mesh vertices
const QgsPoint &vertex = mTriangularMesh.vertices()[vertex_id];
ring.append( vertex );
}
QgsPolygonXY polygon;
polygon.append( ring );
QgsGeometry geom = QgsGeometry::fromPolygonXY( polygon );
QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices() ); //Triangular mesh vertices contains also native mesh vertices
feat.setGeometry( geom );
renderer.renderFeature( feat, mContext );
}

View File

@ -15,6 +15,9 @@
* *
***************************************************************************/
#include <QList>
#include "qgsfeature.h"
#include "qgstriangularmesh.h"
#include "qgsrendercontext.h"
#include "qgscoordinatetransform.h"
@ -63,12 +66,12 @@ static void ENP_centroid( const QPolygonF &pX, double &cx, double &cy )
cy /= ( 6.0 * signedArea );
}
void QgsTriangularMesh::update( QgsMesh *nativeMesh, QgsRenderContext *context )
{
Q_ASSERT( nativeMesh );
Q_ASSERT( context );
mSpatialIndex = QgsSpatialIndex();
mTriangularMesh.vertices.clear();
mTriangularMesh.faces.clear();
mTrianglesToNativeFaces.clear();
@ -143,6 +146,15 @@ void QgsTriangularMesh::update( QgsMesh *nativeMesh, QgsRenderContext *context )
ENP_centroid( poly, cx, cy );
mNativeMeshFaceCentroids[i] = QgsMeshVertex( cx, cy );
}
// CALCULATE SPATIAL INDEX
for ( int i = 0; i < mTriangularMesh.faces.size(); ++i )
{
const QgsMeshFace &face = mTriangularMesh.faces.at( i ) ;
QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
bool success = mSpatialIndex.insertFeature( i, geom.boundingBox() );
Q_UNUSED( success );
}
}
const QVector<QgsMeshVertex> &QgsTriangularMesh::vertices() const
@ -165,3 +177,31 @@ const QVector<int> &QgsTriangularMesh::trianglesToNativeFaces() const
return mTrianglesToNativeFaces;
}
int QgsTriangularMesh::faceIndexForPoint( const QgsPointXY &point ) const
{
const QList<QgsFeatureId> face_indexes = mSpatialIndex.intersects( QgsRectangle( point, point ) );
for ( QgsFeatureId fid : face_indexes )
{
int face_index = static_cast<int>( fid );
const QgsMeshFace &face = mTriangularMesh.faces.at( face_index );
const QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
if ( geom.contains( &point ) )
return face_index;
}
return -1;
}
QgsGeometry QgsMeshUtils::toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
{
QVector<QgsPointXY> ring;
for ( int j = 0; j < face.size(); ++j )
{
int vertex_id = face[j];
Q_ASSERT( vertex_id < vertices.size() );
const QgsPoint &vertex = vertices[vertex_id];
ring.append( vertex );
}
QgsPolygonXY polygon;
polygon.append( ring );
return QgsGeometry::fromPolygonXY( polygon );
}

View File

@ -22,9 +22,10 @@
#define SIP_NO_FILE
#include <QVector>
#include "qgis_core.h"
#include "qgsmeshdataprovider.h"
#include "qgsgeometry.h"
#include "qgsspatialindex.h"
class QgsRenderContext;
@ -40,7 +41,9 @@ struct CORE_EXPORT QgsMesh
/**
* \ingroup core
*
* Triangular/Derived Mesh
* Triangular/Derived Mesh is mesh with vertices in map coordinates. It creates
* spatial index for identification of a triangle that contains a particular point
* on the map.
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
@ -55,7 +58,7 @@ class CORE_EXPORT QgsTriangularMesh
~QgsTriangularMesh() = default;
/**
* Constructs triangular mesh from layer's native mesh and context
* Constructs triangular mesh from layer's native mesh and context. Populates spatial index.
* \param nativeMesh QgsMesh to access native vertices and faces
* \param context Rendering context to estimate number of triagles to create for an face
*/
@ -77,6 +80,12 @@ class CORE_EXPORT QgsTriangularMesh
//! Returns mapping between triangles and original faces
const QVector<int> &trianglesToNativeFaces() const ;
/**
* Returns triangle index that contains the given point, -1 if no such triangle exists
* It uses spatial indexing
*/
int faceIndexForPoint( const QgsPointXY &point ) const ;
private:
// vertices: map CRS; 0-N ... native vertices, N+1 - len ... extra vertices
// faces are derived triangles
@ -85,7 +94,14 @@ class CORE_EXPORT QgsTriangularMesh
// centroids of the native faces in map CRS
QVector<QgsMeshVertex> mNativeMeshFaceCentroids;
QgsSpatialIndex mSpatialIndex;
};
namespace QgsMeshUtils
{
//! Returns face as polygon geometry
QgsGeometry toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices );
};
#endif // QGSTRIANGULARMESH_H