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

qgis · Aug 9, 2018 · c79e1d0 · c79e1d0
1 parent 51b63e6
commit c79e1d0
Show file tree

Hide file tree

Showing 11 changed files with 174 additions and 18 deletions.
diff --git a/python/core/auto_generated/mesh/qgsmeshdataprovider.sip.in b/python/core/auto_generated/mesh/qgsmeshdataprovider.sip.in
@@ -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
 
 };

diff --git a/python/core/auto_generated/mesh/qgsmeshlayer.sip.in b/python/core/auto_generated/mesh/qgsmeshlayer.sip.in
@@ -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

diff --git a/src/core/mesh/qgsmeshdataprovider.cpp b/src/core/mesh/qgsmeshdataprovider.cpp
@@ -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() );

diff --git a/src/core/mesh/qgsmeshdataprovider.h b/src/core/mesh/qgsmeshdataprovider.h
@@ -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;

diff --git a/src/core/mesh/qgsmeshlayer.cpp b/src/core/mesh/qgsmeshlayer.cpp
@@ -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()

diff --git a/src/core/mesh/qgsmeshlayer.h b/src/core/mesh/qgsmeshlayer.h
@@ -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
 
     /**

diff --git a/src/core/mesh/qgsmeshlayerinterpolator.cpp b/src/core/mesh/qgsmeshlayerinterpolator.cpp
@@ -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 ) )

diff --git a/src/core/mesh/qgsmeshlayerinterpolator.h b/src/core/mesh/qgsmeshlayerinterpolator.h
@@ -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;

diff --git a/src/core/mesh/qgsmeshlayerrenderer.cpp b/src/core/mesh/qgsmeshlayerrenderer.cpp
@@ -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 );
   }

diff --git a/src/core/mesh/qgstriangularmesh.cpp b/src/core/mesh/qgstriangularmesh.cpp
@@ -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 );
+}
diff --git a/src/core/mesh/qgstriangularmesh.h b/src/core/mesh/qgstriangularmesh.h
@@ -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