[FEATURE] Rendering of scalar data on mesh layers

Rudimentary support for rendering of scalar data (e.g. water depth)
on mesh map layers.
This commit is contained in:
Martin Dobias 2018-05-10 18:56:08 -04:00 committed by Peter Petrik
parent 3154102aa6
commit 9296528822
12 changed files with 476 additions and 4 deletions

View File

@ -23,6 +23,10 @@ class QgsMeshDatasetValue
QgsMeshDatasetValue is a vector or a scalar value on vertex or face of the mesh with
support of nodata values
.. note::
The API is considered EXPERIMENTAL and can be changed without a notice
.. versionadded:: 3.2
%End
@ -59,6 +63,10 @@ Mesh is a collection of vertices and faces in 2D or 3D space
Base on the underlying data provider/format, whole mesh is either stored in memory or
read on demand
.. note::
The API is considered EXPERIMENTAL and can be changed without a notice
.. versionadded:: 3.2
%End
@ -105,6 +113,10 @@ Dataset is a collection of vector or scalar values on vertices or faces of the
Base on the underlying data provider/format, whole dataset is either stored in memory or
read on demand
.. note::
The API is considered EXPERIMENTAL and can be changed without a notice
.. versionadded:: 3.2
%End
@ -158,6 +170,10 @@ Base class for providing data for :py:class:`QgsMeshLayer`
Responsible for reading native mesh data
.. note::
The API is considered EXPERIMENTAL and can be changed without a notice
.. seealso:: :py:class:`QgsMeshSource`
.. versionadded:: 3.2

View File

@ -64,6 +64,11 @@ is the MDAL connection string. QGIS must be built with MDAL support to allow thi
QString uri = "test/land.2dm";
QgsMeshLayer *scratchLayer = new QgsMeshLayer(uri, "My Scratch Layer", "mdal");
.. note::
The API is considered EXPERIMENTAL and can be changed without a notice
.. versionadded:: 3.2
%End
@ -138,6 +143,8 @@ Toggle rendering of triangular (derived) mesh. Off by default
void setActiveScalarDataset( int index = -1 );
void setActiveVectorDataset( int index = -1 );
int activeScalarDataset() const;
private: // Private methods
QgsMeshLayer( const QgsMeshLayer &rhs );
};

View File

@ -453,6 +453,7 @@ SET(QGIS_CORE_SRCS
mesh/qgsmeshlayerrenderer.cpp
mesh/qgsmeshmemorydataprovider.cpp
mesh/qgstriangularmesh.cpp
mesh/qgsmeshlayerinterpolator.cpp
geometry/qgsabstractgeometry.cpp
geometry/qgsbox3d.cpp
@ -1075,6 +1076,7 @@ SET(QGIS_CORE_HDRS
mesh/qgstriangularmesh.h
mesh/qgsmeshlayerrenderer.h
mesh/qgsmeshlayerinterpolator.h
scalebar/qgsdoubleboxscalebarrenderer.h
scalebar/qgsnumericscalebarrenderer.h

View File

@ -94,6 +94,10 @@ void QgsMeshDatasetValue::setX( double x )
{
mIsNodata = true;
}
else
{
mIsNodata = false;
}
}
void QgsMeshDatasetValue::setY( double y )

View File

@ -43,6 +43,8 @@ typedef QMap<QString, QString> QgsMeshDatasetMetadata;
* QgsMeshDatasetValue is a vector or a scalar value on vertex or face of the mesh with
* support of nodata values
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \since QGIS 3.2
*/
class CORE_EXPORT QgsMeshDatasetValue
@ -83,6 +85,8 @@ class CORE_EXPORT QgsMeshDatasetValue
* Base on the underlying data provider/format, whole mesh is either stored in memory or
* read on demand
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \since QGIS 3.2
*/
class CORE_EXPORT QgsMeshSource SIP_ABSTRACT
@ -123,6 +127,8 @@ class CORE_EXPORT QgsMeshSource SIP_ABSTRACT
* Base on the underlying data provider/format, whole dataset is either stored in memory or
* read on demand
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \since QGIS 3.2
*/
class CORE_EXPORT QgsMeshDatasetSource SIP_ABSTRACT
@ -174,6 +180,8 @@ class CORE_EXPORT QgsMeshDatasetSource SIP_ABSTRACT
*
* Responsible for reading native mesh data
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \see QgsMeshSource
* \since QGIS 3.2
*/

View File

@ -129,6 +129,12 @@ void QgsMeshLayer::toggleTriangularMeshRendering( bool toggle )
void QgsMeshLayer::setActiveScalarDataset( int index )
{
if ( index < 0 )
{
mActiveScalarDataset = -1;
return;
}
Q_ASSERT( dataProvider()->datasetCount() > index );
if ( dataProvider()->datasetHasScalarData( index ) )
mActiveScalarDataset = index;
@ -136,6 +142,12 @@ void QgsMeshLayer::setActiveScalarDataset( int index )
void QgsMeshLayer::setActiveVectorDataset( int index )
{
if ( index < 0 )
{
mActiveVectorDataset = -1;
return;
}
Q_ASSERT( dataProvider()->datasetCount() > index );
if ( !dataProvider()->datasetHasScalarData( index ) )
mActiveVectorDataset = index;

View File

@ -82,6 +82,8 @@ struct QgsMesh;
* QgsMeshLayer *scratchLayer = new QgsMeshLayer(uri, "My Scratch Layer", "mdal");
* \endcode
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \since QGIS 3.2
*/
class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
@ -145,6 +147,8 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
void setActiveScalarDataset( int index = -1 );
void setActiveVectorDataset( int index = -1 );
int activeScalarDataset() const { return mActiveScalarDataset; }
private: // Private methods
/**

View File

@ -0,0 +1,264 @@
/***************************************************************************
qgsmeshlayerinterpolator.cpp
----------------------------
begin : April 2018
copyright : (C) 2018 by Peter Petrik
email : zilolv at gmail dot com
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
///@cond PRIVATE
#include "qgsmeshlayerinterpolator.h"
#include "qgsrasterinterface.h"
#include <QVector2D>
// TODO: use QgsMapToPixel
class MapToPixel
{
public:
MapToPixel( double llX, double llY, double mupp, int rows )
: mLlX( llX ), mLlY( llY ), mMupp( mupp ), mRows( rows ) {}
MapToPixel( const MapToPixel &obj )
: mLlX( obj.mLlX ), mLlY( obj.mLlY ), mMupp( obj.mMupp ), mRows( obj.mRows ) {}
QPointF realToPixel( double rx, double ry ) const
{
double px = ( rx - mLlX ) / mMupp;
double py = mRows - ( ry - mLlY ) / mMupp;
return QPointF( px, py );
}
QPointF realToPixel( const QPointF &pt ) const
{
return realToPixel( pt.x(), pt.y() );
}
QPointF pixelToReal( double px, double py ) const
{
double rx = mLlX + ( px * mMupp );
double ry = mLlY + mMupp * ( mRows - py );
return QPointF( rx, ry );
}
QPointF pixelToReal( const QPointF &pt ) const
{
return pixelToReal( pt.x(), pt.y() );
}
private:
double mLlX, mLlY;
double mMupp; // map units per pixel
double mRows; // (actually integer value)
};
void bbox2rect( const MapToPixel &mtp, const QSize &outputSize, const QgsRectangle &bbox, int &leftLim, int &rightLim, int &topLim, int &bottomLim )
{
QPoint ll = mtp.realToPixel( bbox.xMinimum(), bbox.yMinimum() ).toPoint();
QPoint ur = mtp.realToPixel( bbox.xMaximum(), bbox.yMaximum() ).toPoint();
topLim = std::max( ur.y(), 0 );
bottomLim = std::min( ll.y(), outputSize.height() - 1 );
leftLim = std::max( ll.x(), 0 );
rightLim = std::min( ur.x(), outputSize.width() - 1 );
}
struct MapView
{
MapView(): mtp( 0, 0, 0, 0 ) {}
MapToPixel mtp;
QSize outputSize;
};
static void lam_tol( double &lam )
{
const static double eps = 1e-6;
if ( ( lam < 0.0 ) && ( lam > -eps ) )
{
lam = 0.0;
}
}
bool E3T_physicalToBarycentric( QPointF pA, QPointF pB, QPointF pC, QPointF pP, double &lam1, double &lam2, double &lam3 )
{
if ( pA == pB || pA == pC || pB == pC )
return false; // this is not a valid triangle!
// Compute vectors
QVector2D v0( pC - pA );
QVector2D v1( pB - pA );
QVector2D v2( pP - pA );
// Compute dot products
double dot00 = QVector2D::dotProduct( v0, v0 );
double dot01 = QVector2D::dotProduct( v0, v1 );
double dot02 = QVector2D::dotProduct( v0, v2 );
double dot11 = QVector2D::dotProduct( v1, v1 );
double dot12 = QVector2D::dotProduct( v1, v2 );
// Compute barycentric coordinates
double invDenom = 1.0 / ( dot00 * dot11 - dot01 * dot01 );
lam1 = ( dot11 * dot02 - dot01 * dot12 ) * invDenom;
lam2 = ( dot00 * dot12 - dot01 * dot02 ) * invDenom;
lam3 = 1.0 - lam1 - lam2;
// Apply some tolerance to lam so we can detect correctly border points
lam_tol( lam1 );
lam_tol( lam2 );
lam_tol( lam3 );
// Return if POI is outside triangle
if ( ( lam1 < 0 ) || ( lam2 < 0 ) || ( lam3 < 0 ) )
{
return false;
}
return true;
}
double interpolateFromVerticesData( const QPointF &p1, const QPointF &p2, const QPointF &p3, double val1, double val2, double val3, const QPointF &pt )
{
double lam1, lam2, lam3;
if ( !E3T_physicalToBarycentric( p1, p2, p3, pt, lam1, lam2, lam3 ) )
return std::numeric_limits<double>::quiet_NaN();
return lam1 * val3 + lam2 * val2 + lam3 * val1;
}
double interpolateFromFacesData( const QPointF &p1, const QPointF &p2, const QPointF &p3, double val, const QPointF &pt )
{
double lam1, lam2, lam3;
if ( !E3T_physicalToBarycentric( p1, p2, p3, pt, lam1, lam2, lam3 ) )
return std::numeric_limits<double>::quiet_NaN();
return val;
}
QgsMeshLayerInterpolator::QgsMeshLayerInterpolator( const QgsTriangularMesh &m,
const QVector<double> &datasetValues, bool dataIsOnVertices,
const QgsRenderContext &context )
: mTriangularMesh( m ),
mDatasetValues( datasetValues ),
mContext( context ),
mDataOnVertices( dataIsOnVertices )
{
// figure out image size
QgsRectangle extent = mContext.extent(); // this is extent in layer's coordinate system - but we need it in map coordinate system
QgsMapToPixel mapToPixel = mContext.mapToPixel();
// TODO: what if OTF reprojection is used - see crayfish layer_renderer.py (_calculate_extent)
QgsPointXY topleft = mapToPixel.transform( extent.xMinimum(), extent.yMaximum() );
QgsPointXY bottomright = mapToPixel.transform( extent.xMaximum(), extent.yMinimum() );
int width = bottomright.x() - topleft.x();
int height = bottomright.y() - topleft.y();
mMapView.reset( new MapView() );
mMapView->mtp = MapToPixel( extent.xMinimum(), extent.yMinimum(), mapToPixel.mapUnitsPerPixel(), height );
mMapView->outputSize = QSize( width, height );
}
QgsMeshLayerInterpolator::~QgsMeshLayerInterpolator() = default;
QgsRasterInterface *QgsMeshLayerInterpolator::clone() const
{
return nullptr; // we should not need this (hopefully)
}
Qgis::DataType QgsMeshLayerInterpolator::dataType( int ) const
{
return Qgis::Float32;
}
int QgsMeshLayerInterpolator::bandCount() const
{
return 1;
}
QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
{
QgsRasterBlock *b = new QgsRasterBlock( Qgis::Float32, width, height );
b->setIsNoData(); // assume initially that all values are unset
float *data = reinterpret_cast<float *>( b->bits() );
const QVector<QgsMeshFace> &triangels = mTriangularMesh.triangles();
const QVector<QgsMeshVertex> &vertices = mTriangularMesh.vertices();
// currently expecting that triangulation does not add any new extra vertices on the way
if ( mDataOnVertices )
Q_ASSERT( mDatasetValues.count() == mTriangularMesh.vertices().count() );
for ( int i = 0; i < triangels.size(); ++i )
{
if ( feedback && feedback->isCanceled() )
break;
const QgsMeshFace &face = triangels[i];
const int v1 = face[0], v2 = face[1], v3 = face[2];
const QgsPoint p1 = vertices[v1], p2 = vertices[v2], p3 = vertices[v3];
QgsRectangle bbox;
bbox.combineExtentWith( p1.x(), p1.y() );
bbox.combineExtentWith( p2.x(), p2.y() );
bbox.combineExtentWith( p3.x(), p3.y() );
if ( !extent.intersects( bbox ) )
continue;
// Get the BBox of the element in pixels
int topLim, bottomLim, leftLim, rightLim;
bbox2rect( mMapView->mtp, mMapView->outputSize, bbox, leftLim, rightLim, topLim, bottomLim );
// interpolate in the bounding box of the face
for ( int j = topLim; j <= bottomLim; j++ )
{
float *line = data + ( j * width );
for ( int k = leftLim; k <= rightLim; k++ )
{
double val;
QPointF p = mMapView->mtp.pixelToReal( k, j );
if ( mDataOnVertices )
val = interpolateFromVerticesData(
QPointF( p1.x(), p1.y() ),
QPointF( p2.x(), p2.y() ),
QPointF( p3.x(), p3.y() ),
mDatasetValues[v1],
mDatasetValues[v2],
mDatasetValues[v3],
p );
else
{
int face = mTriangularMesh.trianglesToNativeFaces()[i];
val = interpolateFromFacesData(
QPointF( p1.x(), p1.y() ),
QPointF( p2.x(), p2.y() ),
QPointF( p3.x(), p3.y() ),
mDatasetValues[face],
p
);
}
if ( !qIsNaN( val ) )
{
line[k] = val;
b->setIsData( j, k );
}
}
}
}
return b;
}
///@endcond

View File

@ -0,0 +1,73 @@
/***************************************************************************
qgsmeshlayerinterpolator.h
--------------------------
begin : April 2018
copyright : (C) 2018 by Peter Petrik
email : zilolv at gmail dot com
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#ifndef QGSMESHLAYERINTERPOLATOR_H
#define QGSMESHLAYERINTERPOLATOR_H
class QgsMeshLayer;
class QgsSymbol;
#define SIP_NO_FILE
#include <memory>
#include "qgis.h"
#include "qgsmaplayerrenderer.h"
#include "qgsrendercontext.h"
#include "qgstriangularmesh.h"
#include "qgsrasterinterface.h"
#include "qgssinglebandpseudocolorrenderer.h"
#include "qgsrastershader.h"
#include <QVector2D>
///@cond PRIVATE
struct MapView;
/**
* \ingroup core
* Interpolate mesh scalar dataset to raster block
*
* \since QGIS 3.2
* \note not available in Python bindings
*/
class QgsMeshLayerInterpolator : public QgsRasterInterface
{
public:
QgsMeshLayerInterpolator( const QgsTriangularMesh &m,
const QVector<double> &datasetValues,
bool dataIsOnVertices,
const QgsRenderContext &context );
~QgsMeshLayerInterpolator() override;
virtual QgsRasterInterface *clone() const override;
virtual Qgis::DataType dataType( int ) const override;
virtual int bandCount() const override;
virtual QgsRasterBlock *block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback = nullptr ) override;
private:
const QgsTriangularMesh &mTriangularMesh;
const QVector<double> &mDatasetValues;
const QgsRenderContext &mContext;
std::unique_ptr<MapView> mMapView;
bool mDataOnVertices;
};
///@endcond
#endif // QGSMESHLAYERINTERPOLATOR_H

View File

@ -24,8 +24,9 @@
#include "qgsrenderer.h"
#include "qgssinglesymbolrenderer.h"
#include "qgssymbol.h"
#include "qgssinglebandpseudocolorrenderer.h"
#include "qgsrastershader.h"
#include "qgsmeshlayerinterpolator.h"
QgsMeshLayerRenderer::QgsMeshLayerRenderer( QgsMeshLayer *layer, QgsRenderContext &context )
: QgsMapLayerRenderer( layer->id() )
@ -47,10 +48,46 @@ QgsMeshLayerRenderer::QgsMeshLayerRenderer( QgsMeshLayer *layer, QgsRenderContex
{
mTriangularMeshSymbol.reset( layer->triangularMeshSymbol()->clone() );
}
copyScalarDatasetValues( layer );
}
void QgsMeshLayerRenderer::copyScalarDatasetValues( QgsMeshLayer *layer )
{
int datasetIndex = layer->activeScalarDataset();
if ( datasetIndex != -1 )
{
mDataOnVertices = layer->dataProvider()->datasetIsOnVertices( datasetIndex );
if ( mDataOnVertices )
{
int count = mNativeMesh.vertices.count();
mDatasetValues.resize( count );
for ( int i = 0; i < count; ++i )
{
double v = layer->dataProvider()->datasetValue( datasetIndex, i ).scalar();
mDatasetValues[i] = v;
}
}
else
{
//on faces
int count = mNativeMesh.faces.count();
mDatasetValues.resize( count );
for ( int i = 0; i < count; ++i )
{
double v = layer->dataProvider()->datasetValue( datasetIndex, i ).scalar();
mDatasetValues[i] = v;
}
}
}
}
bool QgsMeshLayerRenderer::render()
{
renderScalarDataset();
renderMesh( mNativeMeshSymbol, mNativeMesh.faces ); // native mesh
renderMesh( mTriangularMeshSymbol, mTriangularMesh.triangles() ); // triangular mesh
@ -91,3 +128,44 @@ void QgsMeshLayerRenderer::renderMesh( const std::unique_ptr<QgsSymbol> &symbol,
renderer.stopRender( mContext );
}
void QgsMeshLayerRenderer::renderScalarDataset()
{
if ( mDatasetValues.isEmpty() )
return;
// figure out image size
QgsRectangle extent = mContext.extent(); // this is extent in layer's coordinate system - but we need it in map coordinate system
QgsMapToPixel mapToPixel = mContext.mapToPixel();
// TODO: what if OTF reprojection is used - see crayfish layer_renderer.py (_calculate_extent)
QgsPointXY topleft = mapToPixel.transform( extent.xMinimum(), extent.yMaximum() );
QgsPointXY bottomright = mapToPixel.transform( extent.xMaximum(), extent.yMinimum() );
int width = bottomright.x() - topleft.x();
int height = bottomright.y() - topleft.y();
double vMin = mDatasetValues[0], vMax = mDatasetValues[0];
for ( int i = 1; i < mDatasetValues.count(); ++i )
{
double v = mDatasetValues[i];
if ( v < vMin ) vMin = v;
if ( v > vMax ) vMax = v;
}
QList<QgsColorRampShader::ColorRampItem> lst;
lst << QgsColorRampShader::ColorRampItem( vMin, Qt::red, QString::number( vMin ) );
lst << QgsColorRampShader::ColorRampItem( vMax, Qt::blue, QString::number( vMax ) );
QgsColorRampShader *fcn = new QgsColorRampShader( vMin, vMax );
fcn->setColorRampItemList( lst );
QgsRasterShader *sh = new QgsRasterShader( 0, 1000 );
sh->setRasterShaderFunction( fcn ); // takes ownership of fcn
QgsMeshLayerInterpolator interpolator( mTriangularMesh, mDatasetValues, mDataOnVertices, mContext );
QgsSingleBandPseudoColorRenderer r( &interpolator, 0, sh ); // takes ownership of sh
QgsRasterBlock *bl = r.block( 0, extent, width, height ); // TODO: feedback
Q_ASSERT( bl );
QImage img = bl->image();
mContext.painter()->drawImage( 0, 0, img );
delete bl;
}

View File

@ -49,7 +49,8 @@ class QgsMeshLayerRenderer : public QgsMapLayerRenderer
private:
void renderMesh( const std::unique_ptr<QgsSymbol> &symbol, const QVector<QgsMeshFace> &faces );
void renderScalarDataset();
void copyScalarDatasetValues( QgsMeshLayer *layer );
protected:
// copy from mesh layer
@ -59,7 +60,8 @@ class QgsMeshLayerRenderer : public QgsMapLayerRenderer
QgsTriangularMesh mTriangularMesh;
// copy of the scalar dataset
QVector<double> mDatasetValues;
bool mDataOnVertices;
// copy from mesh layer
std::unique_ptr<QgsSymbol> mNativeMeshSymbol = nullptr;

View File

@ -42,6 +42,8 @@ struct CORE_EXPORT QgsMesh
*
* Triangular/Derived Mesh
*
* \note The API is considered EXPERIMENTAL and can be changed without a notice
*
* \since QGIS 3.2
*/
class CORE_EXPORT QgsTriangularMesh