Skip to content

Commit 0d9a15b

Browse files
vcloarecvcloarec
authoredMay 5, 2020
Snap on mesh elements (#36171)
[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 :
·
ltr-3_40final-3_14_0
1 parent c22c73d commit 0d9a15b

File tree

4 files changed

+260
-45
lines changed

4 files changed

+260
-45
lines changed
 

‎python/core/auto_generated/mesh/qgsmeshlayer.sip.in

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -317,6 +317,34 @@ Sets the reference time of the layer
317317

318318
:param referenceTime: the reference time
319319

320+
.. versionadded:: 3.14
321+
%End
322+
323+
QgsPointXY snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius );
324+
%Docstring
325+
Returns the position of the snapped point on the mesh element closest to ``point`` intersecting with
326+
the searching area defined by ``point`` and ``searchRadius``
327+
328+
For vertex, the snapped position is the vertex position
329+
For edge, the snapped position is the projected point on the edge, extremity of edge if outside the edge
330+
For face, the snapped position is the centroid of the face
331+
The returned position is in map coordinates.
332+
333+
.. note::
334+
335+
It uses previously cached and indexed triangular mesh
336+
and so if the layer has not been rendered previously
337+
(e.g. when used in a script) it returns empty :py:class:`QgsPointXY`
338+
339+
.. seealso:: :py:func:`updateTriangularMesh`
340+
341+
:param elementType: the type of element to snap
342+
:param point: the center of the search area in map coordinates
343+
:param searchRadius: the radius of the search area in map units
344+
345+
:return: the position of the snapped point on the closest element, empty QgsPointXY if no element of type ``elementType``
346+
347+
320348
.. versionadded:: 3.14
321349
%End
322350

‎src/core/mesh/qgsmeshlayer.cpp

Lines changed: 152 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -354,62 +354,38 @@ QgsMesh3dDataBlock QgsMeshLayer::dataset3dValue( const QgsMeshDatasetIndex &inde
354354
QgsMeshDatasetValue QgsMeshLayer::dataset1dValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius ) const
355355
{
356356
QgsMeshDatasetValue value;
357-
QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
357+
QgsPointXY projectedPoint;
358+
int selectedIndex = closestEdge( point, searchRadius, projectedPoint );
358359
const QgsTriangularMesh *mesh = triangularMesh();
359-
if ( mesh &&
360-
mesh->contains( QgsMesh::Edge ) &&
361-
mDataProvider->isValid() &&
362-
index.isValid() )
360+
if ( selectedIndex >= 0 )
363361
{
364-
// search for the closest edge in rectangle from point
365-
const QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( searchRectangle );
366-
int selectedIndex = -1;
367-
double sqrMaxDistFromPoint = pow( searchRadius, 2 );
368-
QgsPointXY projectedPoint;
369-
for ( const int edgeIndex : edgeIndexes )
362+
const QgsMeshDatasetGroupMetadata::DataType dataType = dataProvider()->datasetGroupMetadata( index ).dataType();
363+
switch ( dataType )
370364
{
371-
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
372-
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
373-
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
374-
QgsPointXY projPoint;
375-
double sqrDist = point.sqrDistToSegment( vertex1.x(), vertex1.y(), vertex2.x(), vertex2.y(), projPoint );
376-
if ( sqrDist < sqrMaxDistFromPoint )
365+
case QgsMeshDatasetGroupMetadata::DataOnEdges:
377366
{
378-
selectedIndex = edgeIndex;
379-
projectedPoint = projPoint;
380-
sqrMaxDistFromPoint = sqrDist;
367+
value = dataProvider()->datasetValue( index, selectedIndex );
381368
}
382-
}
369+
break;
383370

384-
if ( selectedIndex >= 0 )
385-
{
386-
const QgsMeshDatasetGroupMetadata::DataType dataType = dataProvider()->datasetGroupMetadata( index ).dataType();
387-
switch ( dataType )
371+
case QgsMeshDatasetGroupMetadata::DataOnVertices:
388372
{
389-
case QgsMeshDatasetGroupMetadata::DataOnEdges:
390-
{
391-
value = dataProvider()->datasetValue( index, selectedIndex );
392-
}
393-
break;
394-
395-
case QgsMeshDatasetGroupMetadata::DataOnVertices:
396-
{
397-
const QgsMeshEdge &edge = mesh->edges()[selectedIndex];
398-
const int v1 = edge.first, v2 = edge.second;
399-
const QgsPoint p1 = mesh->vertices()[v1], p2 = mesh->vertices()[v2];
400-
const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
401-
const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
402-
double edgeLength = p1.distance( p2 );
403-
double dist1 = p1.distance( projectedPoint.x(), projectedPoint.y() );
404-
value = QgsMeshLayerUtils::interpolateFromVerticesData( dist1 / edgeLength, val1, val2 );
405-
}
406-
break;
407-
default:
408-
break;
373+
const QgsMeshEdge &edge = mesh->edges()[selectedIndex];
374+
const int v1 = edge.first, v2 = edge.second;
375+
const QgsPoint p1 = mesh->vertices()[v1], p2 = mesh->vertices()[v2];
376+
const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
377+
const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
378+
double edgeLength = p1.distance( p2 );
379+
double dist1 = p1.distance( projectedPoint.x(), projectedPoint.y() );
380+
value = QgsMeshLayerUtils::interpolateFromVerticesData( dist1 / edgeLength, val1, val2 );
409381
}
382+
break;
383+
default:
384+
break;
410385
}
411386
}
412387

388+
413389
return value;
414390
}
415391

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

472448
}
473449

450+
int QgsMeshLayer::closestEdge( const QgsPointXY &point, double searchRadius, QgsPointXY &projectedPoint ) const
451+
{
452+
QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
453+
const QgsTriangularMesh *mesh = triangularMesh();
454+
// search for the closest edge in search area from point
455+
const QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( searchRectangle );
456+
int selectedIndex = -1;
457+
if ( mesh &&
458+
mesh->contains( QgsMesh::Edge ) &&
459+
mDataProvider->isValid() )
460+
{
461+
double sqrMaxDistFromPoint = pow( searchRadius, 2 );
462+
for ( const int edgeIndex : edgeIndexes )
463+
{
464+
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
465+
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
466+
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
467+
QgsPointXY projPoint;
468+
double sqrDist = point.sqrDistToSegment( vertex1.x(), vertex1.y(), vertex2.x(), vertex2.y(), projPoint );
469+
if ( sqrDist < sqrMaxDistFromPoint )
470+
{
471+
selectedIndex = edgeIndex;
472+
projectedPoint = projPoint;
473+
sqrMaxDistFromPoint = sqrDist;
474+
}
475+
}
476+
}
477+
478+
return selectedIndex;
479+
}
480+
474481
QgsMeshDatasetIndex QgsMeshLayer::staticVectorDatasetIndex() const
475482
{
476483
return mStaticVectorDatasetIndex;
@@ -484,6 +491,106 @@ void QgsMeshLayer::setReferenceTime( const QDateTime &referenceTime )
484491
temporalProperties()->setReferenceTime( referenceTime, nullptr );
485492
}
486493

494+
QgsPointXY QgsMeshLayer::snapOnVertex( const QgsPointXY &point, double searchRadius )
495+
{
496+
const QgsTriangularMesh *mesh = triangularMesh();
497+
QgsPointXY exactPosition;
498+
if ( !mesh )
499+
return exactPosition;
500+
QgsRectangle rectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
501+
double maxDistance = searchRadius;
502+
//attempt to snap on edges's vertices
503+
QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( rectangle );
504+
for ( const int edgeIndex : edgeIndexes )
505+
{
506+
const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
507+
const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
508+
const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
509+
double dist1 = point.distance( vertex1 );
510+
double dist2 = point.distance( vertex2 );
511+
if ( dist1 < maxDistance )
512+
{
513+
maxDistance = dist1;
514+
exactPosition = vertex1;
515+
}
516+
if ( dist2 < maxDistance )
517+
{
518+
maxDistance = dist2;
519+
exactPosition = vertex2;
520+
}
521+
}
522+
523+
//attempt to snap on face's vertices
524+
QList<int> faceIndexes = mesh->faceIndexesForRectangle( rectangle );
525+
for ( const int faceIndex : faceIndexes )
526+
{
527+
const QgsMeshFace &face = mesh->triangles().at( faceIndex );
528+
for ( int i = 0; i < 3; ++i )
529+
{
530+
const QgsMeshVertex &vertex = mesh->vertices()[face.at( i )];
531+
double dist = point.distance( vertex );
532+
if ( dist < maxDistance )
533+
{
534+
maxDistance = dist;
535+
exactPosition = vertex;
536+
}
537+
}
538+
}
539+
540+
return exactPosition;
541+
}
542+
543+
QgsPointXY QgsMeshLayer::snapOnEdge( const QgsPointXY &point, double searchRadius )
544+
{
545+
QgsPointXY projectedPoint;
546+
closestEdge( point, searchRadius, projectedPoint );
547+
548+
return projectedPoint;
549+
}
550+
551+
QgsPointXY QgsMeshLayer::snapOnFace( const QgsPointXY &point, double searchRadius )
552+
{
553+
const QgsTriangularMesh *mesh = triangularMesh();
554+
QgsPointXY centroidPosition;
555+
if ( !mesh )
556+
return centroidPosition;
557+
QgsRectangle rectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
558+
double maxDistance = std::numeric_limits<double>::max();
559+
560+
QList<int> faceIndexes = mesh->faceIndexesForRectangle( rectangle );
561+
for ( const int faceIndex : faceIndexes )
562+
{
563+
int nativefaceIndex = mesh->trianglesToNativeFaces().at( faceIndex );
564+
if ( nativefaceIndex < 0 && nativefaceIndex >= mesh->faceCentroids().count() )
565+
continue;
566+
const QgsPointXY centroid = mesh->faceCentroids()[nativefaceIndex];
567+
double dist = point.distance( centroid );
568+
if ( dist < maxDistance )
569+
{
570+
maxDistance = dist;
571+
centroidPosition = centroid;
572+
}
573+
}
574+
575+
return centroidPosition;
576+
}
577+
578+
QgsPointXY QgsMeshLayer::snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius )
579+
{
580+
switch ( elementType )
581+
{
582+
case QgsMesh::Vertex:
583+
return snapOnVertex( point, searchRadius );
584+
break;
585+
case QgsMesh::Edge:
586+
return snapOnEdge( point, searchRadius );
587+
break;
588+
case QgsMesh::Face:
589+
return snapOnFace( point, searchRadius );
590+
break;
591+
}
592+
}
593+
487594
QgsMeshDatasetIndex QgsMeshLayer::staticScalarDatasetIndex() const
488595
{
489596
return mStaticScalarDatasetIndex;

‎src/core/mesh/qgsmeshlayer.h

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -377,6 +377,29 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
377377
*/
378378
void setReferenceTime( const QDateTime &referenceTime );
379379

380+
/**
381+
* Returns the position of the snapped point on the mesh element closest to \a point intersecting with
382+
* the searching area defined by \a point and \a searchRadius
383+
*
384+
* For vertex, the snapped position is the vertex position
385+
* For edge, the snapped position is the projected point on the edge, extremity of edge if outside the edge
386+
* For face, the snapped position is the centroid of the face
387+
* The returned position is in map coordinates.
388+
*
389+
* \note It uses previously cached and indexed triangular mesh
390+
* and so if the layer has not been rendered previously
391+
* (e.g. when used in a script) it returns empty QgsPointXY
392+
* \see updateTriangularMesh
393+
*
394+
* \param elementType the type of element to snap
395+
* \param point the center of the search area in map coordinates
396+
* \param searchRadius the radius of the search area in map units
397+
* \return the position of the snapped point on the closest element, empty QgsPointXY if no element of type \a elementType
398+
*
399+
* \since QGIS 3.14
400+
*/
401+
QgsPointXY snapOnElement( QgsMesh::ElementType elementType, const QgsPointXY &point, double searchRadius );
402+
380403
public slots:
381404

382405
/**
@@ -467,6 +490,17 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
467490

468491
QgsMeshDatasetIndex mStaticScalarDatasetIndex;
469492
QgsMeshDatasetIndex mStaticVectorDatasetIndex;
493+
494+
int closestEdge( const QgsPointXY &point, double searchRadius, QgsPointXY &projectedPoint ) const;
495+
496+
//! Returns the exact position in map coordinates of the closest vertex in the search area
497+
QgsPointXY snapOnVertex( const QgsPointXY &point, double searchRadius );
498+
499+
//!Returns the position of the projected point on the closest edge in the search area
500+
QgsPointXY snapOnEdge( const QgsPointXY &point, double searchRadius );
501+
502+
//!Returns the position of the centroid point on the closest face in the search area
503+
QgsPointXY snapOnFace( const QgsPointXY &point, double searchRadius );
470504
};
471505

472506
#endif //QGSMESHLAYER_H

‎tests/src/core/testqgsmeshlayer.cpp

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,7 @@ class TestQgsMeshLayer : public QObject
8383

8484
void test_mesh_simplification();
8585

86+
void test_snap_on_mesh();
8687
void test_dataset_value_from_layer();
8788
};
8889

@@ -911,6 +912,51 @@ void TestQgsMeshLayer::test_mesh_simplification()
911912
delete m;
912913
}
913914

915+
void TestQgsMeshLayer::test_snap_on_mesh()
916+
{
917+
//1D mesh
918+
mMdal1DLayer->updateTriangularMesh();
919+
double searchRadius = 10;
920+
921+
QgsPointXY snappedPoint;
922+
923+
//1D mesh
924+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
925+
QCOMPARE( snappedPoint, QgsPointXY() );
926+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 1002, 2005 ), searchRadius );
927+
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
928+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 1002, 2005 ), searchRadius );
929+
QCOMPARE( snappedPoint, QgsPointXY( 1002, 2000 ) );
930+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 998, 2005 ), searchRadius );
931+
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
932+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Edge, QgsPointXY( 990, 2010 ), searchRadius );
933+
QCOMPARE( snappedPoint, QgsPointXY() );
934+
snappedPoint = mMdal1DLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 2002, 2998 ), searchRadius );
935+
QCOMPARE( snappedPoint, QgsPointXY( 2000, 3000 ) );
936+
937+
938+
//2D mesh
939+
mMdalLayer->updateTriangularMesh();
940+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
941+
QCOMPARE( snappedPoint, QgsPointXY() );
942+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY(), searchRadius );
943+
QCOMPARE( snappedPoint, QgsPointXY() );
944+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 1002, 2005 ), searchRadius );
945+
QCOMPARE( snappedPoint, QgsPointXY( 1000, 2000 ) );
946+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Vertex, QgsPointXY( 2002, 2998 ), searchRadius );
947+
QCOMPARE( snappedPoint, QgsPointXY( 2000, 3000 ) );
948+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 998, 1998 ), searchRadius );
949+
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
950+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 1002, 2001 ), searchRadius );
951+
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
952+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 1998, 2998 ), searchRadius );
953+
QCOMPARE( snappedPoint, QgsPointXY( 1500, 2500 ) );
954+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 2002, 1998 ), searchRadius );
955+
QCOMPARE( snappedPoint, QgsPointXY( 2333.33333333, 2333.333333333 ) );
956+
snappedPoint = mMdalLayer->snapOnElement( QgsMesh::Face, QgsPointXY( 500, 500 ), searchRadius );
957+
QCOMPARE( snappedPoint, QgsPointXY() );
958+
}
959+
914960
void TestQgsMeshLayer::test_dataset_value_from_layer()
915961
{
916962
QgsMeshDatasetValue value;

0 commit comments

Comments
 (0)
Please sign in to comment.