Skip to content

Commit

Permalink
Snap on mesh elements (#36171)
Browse files Browse the repository at this point in the history 
[API] [FEATURE] Add a method in API to snap on mesh elements (vertex, edge, face)
The method returns the position of the snapped point on the closest mesh element :
  • Loading branch information
@vcloarec
vcloarec committed May 5, 2020
1 parent c22c73d commit 0d9a15b
Show file tree
Hide file tree
Showing 4 changed files with 260 additions and 45 deletions.
28 changes: 28 additions & 0 deletions python/core/auto_generated/mesh/qgsmeshlayer.sip.in
View file
Expand Up @@ -317,6 +317,34 @@ Sets the reference time of the layer

:param referenceTime: the reference time

.. versionadded:: 3.14
%End

QgsPointXY snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius );
%Docstring
Returns the position of the snapped point on the mesh element closest to ``point`` intersecting with
the searching area defined by ``point`` and ``searchRadius``

For vertex, the snapped position is the vertex position
For edge, the snapped position is the projected point on the edge, extremity of edge if outside the edge
For face, the snapped position is the centroid of the face
The returned position is in map coordinates.

.. note::

It uses previously cached and indexed triangular mesh
and so if the layer has not been rendered previously
(e.g. when used in a script) it returns empty :py:class:`QgsPointXY`

.. seealso:: :py:func:`updateTriangularMesh`

:param elementType: the type of element to snap
:param point: the center of the search area in map coordinates
:param searchRadius: the radius of the search area in map units

:return: the position of the snapped point on the closest element, empty QgsPointXY if no element of type ``elementType``


.. versionadded:: 3.14
%End

Expand Down
197 changes: 152 additions & 45 deletions src/core/mesh/qgsmeshlayer.cpp
View file
Expand Up @@ -354,62 +354,38 @@ QgsMesh3dDataBlock QgsMeshLayer::dataset3dValue( const QgsMeshDatasetIndex &inde
QgsMeshDatasetValue QgsMeshLayer::dataset1dValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius ) const
{
QgsMeshDatasetValue value;
QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
QgsPointXY projectedPoint;
int selectedIndex = closestEdge( point, searchRadius, projectedPoint );
const QgsTriangularMesh *mesh = triangularMesh();
if ( mesh &&
mesh->contains( QgsMesh::Edge ) &&
mDataProvider->isValid() &&
index.isValid() )
if ( selectedIndex >= 0 )
{
// search for the closest edge in rectangle from point
const QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( searchRectangle );
int selectedIndex = -1;
double sqrMaxDistFromPoint = pow( searchRadius, 2 );
QgsPointXY projectedPoint;
for ( const int edgeIndex : edgeIndexes )
const QgsMeshDatasetGroupMetadata::DataType dataType = dataProvider()->datasetGroupMetadata( index ).dataType();
switch ( dataType )
{
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
QgsPointXY projPoint;
double sqrDist = point.sqrDistToSegment( vertex1.x(), vertex1.y(), vertex2.x(), vertex2.y(), projPoint );
if ( sqrDist < sqrMaxDistFromPoint )
case QgsMeshDatasetGroupMetadata::DataOnEdges:
{
selectedIndex = edgeIndex;
projectedPoint = projPoint;
sqrMaxDistFromPoint = sqrDist;
value = dataProvider()->datasetValue( index, selectedIndex );
}
}
break;

if ( selectedIndex >= 0 )
{
const QgsMeshDatasetGroupMetadata::DataType dataType = dataProvider()->datasetGroupMetadata( index ).dataType();
switch ( dataType )
case QgsMeshDatasetGroupMetadata::DataOnVertices:
{
case QgsMeshDatasetGroupMetadata::DataOnEdges:
{
value = dataProvider()->datasetValue( index, selectedIndex );
}
break;

case QgsMeshDatasetGroupMetadata::DataOnVertices:
{
const QgsMeshEdge &edge = mesh->edges()[selectedIndex];
const int v1 = edge.first, v2 = edge.second;
const QgsPoint p1 = mesh->vertices()[v1], p2 = mesh->vertices()[v2];
const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
double edgeLength = p1.distance( p2 );
double dist1 = p1.distance( projectedPoint.x(), projectedPoint.y() );
value = QgsMeshLayerUtils::interpolateFromVerticesData( dist1 / edgeLength, val1, val2 );
}
break;
default:
break;
const QgsMeshEdge &edge = mesh->edges()[selectedIndex];
const int v1 = edge.first, v2 = edge.second;
const QgsPoint p1 = mesh->vertices()[v1], p2 = mesh->vertices()[v2];
const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
double edgeLength = p1.distance( p2 );
double dist1 = p1.distance( projectedPoint.x(), projectedPoint.y() );
value = QgsMeshLayerUtils::interpolateFromVerticesData( dist1 / edgeLength, val1, val2 );
}
break;
default:
break;
}
}


return value;
}

Expand Down Expand Up @@ -471,6 +447,37 @@ void QgsMeshLayer::onDatasetGroupsAdded( int count )

}

int QgsMeshLayer::closestEdge( const QgsPointXY &point, double searchRadius, QgsPointXY &projectedPoint ) const
{
QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
const QgsTriangularMesh *mesh = triangularMesh();
// search for the closest edge in search area from point
const QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( searchRectangle );
int selectedIndex = -1;
if ( mesh &&
mesh->contains( QgsMesh::Edge ) &&
mDataProvider->isValid() )
{
double sqrMaxDistFromPoint = pow( searchRadius, 2 );
for ( const int edgeIndex : edgeIndexes )
{
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
QgsPointXY projPoint;
double sqrDist = point.sqrDistToSegment( vertex1.x(), vertex1.y(), vertex2.x(), vertex2.y(), projPoint );
if ( sqrDist < sqrMaxDistFromPoint )
{
selectedIndex = edgeIndex;
projectedPoint = projPoint;
sqrMaxDistFromPoint = sqrDist;
}
}
}

return selectedIndex;
}

QgsMeshDatasetIndex QgsMeshLayer::staticVectorDatasetIndex() const
{
return mStaticVectorDatasetIndex;
Expand All @@ -484,6 +491,106 @@ void QgsMeshLayer::setReferenceTime( const QDateTime &referenceTime )
temporalProperties()->setReferenceTime( referenceTime, nullptr );
}

QgsPointXY QgsMeshLayer::snapOnVertex( const QgsPointXY &point, double searchRadius )
{
const QgsTriangularMesh *mesh = triangularMesh();
QgsPointXY exactPosition;
if ( !mesh )
return exactPosition;
QgsRectangle rectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
double maxDistance = searchRadius;
//attempt to snap on edges's vertices
QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( rectangle );
for ( const int edgeIndex : edgeIndexes )
{
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
double dist1 = point.distance( vertex1 );
double dist2 = point.distance( vertex2 );
if ( dist1 < maxDistance )
{
maxDistance = dist1;
exactPosition = vertex1;
}
if ( dist2 < maxDistance )
{
maxDistance = dist2;
exactPosition = vertex2;
}
}

//attempt to snap on face's vertices
QList<int> faceIndexes = mesh->faceIndexesForRectangle( rectangle );
for ( const int faceIndex : faceIndexes )
{
const QgsMeshFace &face = mesh->triangles().at( faceIndex );
for ( int i = 0; i < 3; ++i )
{
const QgsMeshVertex &vertex = mesh->vertices()[face.at( i )];
double dist = point.distance( vertex );
if ( dist < maxDistance )
{
maxDistance = dist;
exactPosition = vertex;
}
}
}

return exactPosition;
}

QgsPointXY QgsMeshLayer::snapOnEdge( const QgsPointXY &point, double searchRadius )
{
QgsPointXY projectedPoint;
closestEdge( point, searchRadius, projectedPoint );

return projectedPoint;
}

QgsPointXY QgsMeshLayer::snapOnFace( const QgsPointXY &point, double searchRadius )
{
const QgsTriangularMesh *mesh = triangularMesh();
QgsPointXY centroidPosition;
if ( !mesh )
return centroidPosition;
QgsRectangle rectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
double maxDistance = std::numeric_limits<double>::max();

QList<int> faceIndexes = mesh->faceIndexesForRectangle( rectangle );
for ( const int faceIndex : faceIndexes )
{
int nativefaceIndex = mesh->trianglesToNativeFaces().at( faceIndex );
if ( nativefaceIndex < 0 && nativefaceIndex >= mesh->faceCentroids().count() )
continue;
const QgsPointXY centroid = mesh->faceCentroids()[nativefaceIndex];
double dist = point.distance( centroid );
if ( dist < maxDistance )
{
maxDistance = dist;
centroidPosition = centroid;
}
}

return centroidPosition;
}

QgsPointXY QgsMeshLayer::snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius )
{
switch ( elementType )
{
case QgsMesh::Vertex:
return snapOnVertex( point, searchRadius );
break;
case QgsMesh::Edge:
return snapOnEdge( point, searchRadius );
break;
case QgsMesh::Face:
return snapOnFace( point, searchRadius );
break;
}
}

QgsMeshDatasetIndex QgsMeshLayer::staticScalarDatasetIndex() const
{
return mStaticScalarDatasetIndex;
Expand Down
34 changes: 34 additions & 0 deletions src/core/mesh/qgsmeshlayer.h
View file
Expand Up @@ -377,6 +377,29 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
*/
void setReferenceTime( const QDateTime &referenceTime );

/**
* Returns the position of the snapped point on the mesh element closest to \a point intersecting with
* the searching area defined by \a point and \a searchRadius
*
* For vertex, the snapped position is the vertex position
* For edge, the snapped position is the projected point on the edge, extremity of edge if outside the edge
* For face, the snapped position is the centroid of the face
* The returned position is in map coordinates.
*
* \note It uses previously cached and indexed triangular mesh
* and so if the layer has not been rendered previously
* (e.g. when used in a script) it returns empty QgsPointXY
* \see updateTriangularMesh
*
* \param elementType the type of element to snap
* \param point the center of the search area in map coordinates
* \param searchRadius the radius of the search area in map units
* \return the position of the snapped point on the closest element, empty QgsPointXY if no element of type \a elementType
*
* \since QGIS 3.14
*/
QgsPointXY snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius );

public slots:

/**
Expand Down Expand Up @@ -467,6 +490,17 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer

QgsMeshDatasetIndex mStaticScalarDatasetIndex;
QgsMeshDatasetIndex mStaticVectorDatasetIndex;

int closestEdge( const QgsPointXY &point, double searchRadius, QgsPointXY &projectedPoint ) const;

//! Returns the exact position in map coordinates of the closest vertex in the search area
QgsPointXY snapOnVertex( const QgsPointXY &point, double searchRadius );

//!Returns the position of the projected point on the closest edge in the search area
QgsPointXY snapOnEdge( const QgsPointXY &point, double searchRadius );

//!Returns the position of the centroid point on the closest face in the search area
QgsPointXY snapOnFace( const QgsPointXY &point, double searchRadius );
};

#endif //QGSMESHLAYER_H
46 changes: 46 additions & 0 deletions tests/src/core/testqgsmeshlayer.cpp
View file
Expand Up @@ -83,6 +83,7 @@ class TestQgsMeshLayer : public QObject

void test_mesh_simplification();

void test_snap_on_mesh();
void test_dataset_value_from_layer();
};

Expand Down Expand Up @@ -911,6 +912,51 @@ void TestQgsMeshLayer::test_mesh_simplification()
delete m;
}

void TestQgsMeshLayer::test_snap_on_mesh()
{
//1D mesh
mMdal1DLayer->updateTriangularMesh();
double searchRadius = 10;

QgsPointXY snappedPoint;

//1D mesh
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY() );
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 1002, 2005 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 1002, 2005 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1002, 2000 ) );
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 998, 2005 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 990, 2010 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY() );
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 2002, 2998 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 2000, 3000 ) );


//2D mesh
mMdalLayer->updateTriangularMesh();
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY() );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY() );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 1002, 2005 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 2002, 2998 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 2000, 3000 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 998, 1998 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 1002, 2001 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 1998, 2998 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 2002, 1998 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY( 2333.33333333, 2333.333333333 ) );
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 500, 500 ), searchRadius );
QCOMPARE( snappedPoint, QgsPointXY() );
}

void TestQgsMeshLayer::test_dataset_value_from_layer()
{
QgsMeshDatasetValue value;
Expand Down

0 comments on commit 0d9a15b

Please sign in to comment.