Skip to content
Permalink
Browse files

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

  • Loading branch information
PeterPetrik committed Aug 9, 2018
1 parent 51b63e6 commit c79e1d0601d6b79358e02f05893d1b8689c30f86
@@ -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

};
@@ -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
@@ -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() );
@@ -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;
@@ -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()
@@ -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

/**
@@ -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 ) )
@@ -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;
@@ -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 );
}
@@ -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 );
}
@@ -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

0 comments on commit c79e1d0

Please sign in to comment.
You can’t perform that action at this time.