Mesh editing delaunay refinement

python/PyQt6/core/auto_generated/mesh/qgsmesheditor.sip.in

@@ -260,6 +260,14 @@ Returns the count of valid vertices, that is non void vertices in the mesh
     int maximumVerticesPerFace() const;
 %Docstring
 Returns the maximum count of vertices per face that the mesh can support
+%End
+
+    void addVertexWithDelaunayRefinement( const QgsMeshVertex &vertex, const double tolerance );
+%Docstring
+Add a vertex in a face with Delaunay refinement of neighboring faces
+All neighboring faces sharing a vertex will be refined to satisfy the Delaunay condition
+
+.. versionadded:: 3.42
 %End
 
   signals:

python/core/auto_generated/mesh/qgsmesheditor.sip.in

@@ -260,6 +260,14 @@ Returns the count of valid vertices, that is non void vertices in the mesh
     int maximumVerticesPerFace() const;
 %Docstring
 Returns the maximum count of vertices per face that the mesh can support
+%End
+
+    void addVertexWithDelaunayRefinement( const QgsMeshVertex &vertex, const double tolerance );
+%Docstring
+Add a vertex in a face with Delaunay refinement of neighboring faces
+All neighboring faces sharing a vertex will be refined to satisfy the Delaunay condition
+
+.. versionadded:: 3.42
 %End
 
   signals:

src/app/mesh/qgsmaptooleditmeshframe.cpp

@@ -238,20 +238,29 @@ QgsMeshEditDigitizingAction::QgsMeshEditDigitizingAction( QObject *parent )
   int interpolateFromValue = settings.enumValue( QStringLiteral( "UI/Mesh/zValueFrom" ), PreferMeshThenZWidget );
   mComboZValueType->setCurrentIndex( interpolateFromValue );
 
-  gLayout->addWidget( labelZValueType, 2, 0, 1, 3 );
-  gLayout->addWidget( mComboZValueType, 2, 3, 1, 1 );
+  mCheckBoxRefineNeighboringFaces = new QCheckBox( tr( "Refine neighboring faces when adding vertices" ) );
+  mCheckBoxRefineNeighboringFaces->setToolTip( "Flip edges that do not fulfil delaunay rule on triangular faces that share at least one vertex with the face that new vertex was added to." );
 
+  bool refineNeighboringFaces = settings.value( QStringLiteral( "UI/Mesh/refineNeighboringFaces" ) ).toBool();
+  mCheckBoxRefineNeighboringFaces->setChecked( refineNeighboringFaces );
+
+  gLayout->addWidget( labelZValueType, 1, 0, 1, 1 );
+  gLayout->addWidget( mComboZValueType, 1, 1, 1, 1 );
+  gLayout->addWidget( mCheckBoxRefineNeighboringFaces, 2, 0, 1, 2 );
 
   QWidget *w = new QWidget();
   w->setLayout( gLayout );
   setDefaultWidget( w );
+
+  connect( mCheckBoxRefineNeighboringFaces, &QCheckBox::toggled, this, &QgsMeshEditDigitizingAction::updateSettings );
 }
 
 void QgsMeshEditDigitizingAction::updateSettings()
 {
   QgsSettings settings;
 
   settings.setEnumValue( QStringLiteral( "UI/Mesh/zValueFrom" ), static_cast<ZValueSource>( mComboZValueType->currentData().toInt() ) );
+  settings.setValue( QStringLiteral( "UI/Mesh/refineNeighboringFaces" ), mCheckBoxRefineNeighboringFaces->isChecked() );
 }
 
 QgsMeshEditDigitizingAction::ZValueSource QgsMeshEditDigitizingAction::zValueSourceType() const
@@ -264,6 +273,11 @@ void QgsMeshEditDigitizingAction::setZValueType( QgsMeshEditDigitizingAction::ZV
   mComboZValueType->setCurrentIndex( mComboZValueType->findData( zValueSource ) );
 }
 
+bool QgsMeshEditDigitizingAction::refineNeighboringFaces() const
+{
+  return mCheckBoxRefineNeighboringFaces->isChecked();
+}
+
 //
 // QgsMapToolEditMeshFrame
 //
@@ -2798,10 +2812,18 @@ void QgsMapToolEditMeshFrame::addVertex(
   }
 
   const QVector<QgsMeshVertex> points( 1, QgsMeshVertex( effectivePoint.x(), effectivePoint.y(), zValue ) );
+
   if ( mCurrentEditor )
   {
     double tolerance = QgsTolerance::vertexSearchRadius( canvas()->mapSettings() );
-    mCurrentEditor->addVertices( points, tolerance );
+    if ( mWidgetActionDigitizing->refineNeighboringFaces() && mCurrentFaceIndex != -1 )
+    {
+      mCurrentEditor->addVertexWithDelaunayRefinement( points.first(), tolerance );
+    }
+    else
+    {
+      mCurrentEditor->addVertices( points, tolerance );
+    }
   }
 }
 

src/app/mesh/qgsmaptooleditmeshframe.h

@@ -141,13 +141,19 @@ class QgsMeshEditDigitizingAction : public QWidgetAction
     //! Returns type of z value obtaining
     QgsMeshEditDigitizingAction::ZValueSource zValueSourceType() const;
 
+    //! Returns if neighboring faces should be refined when adding vertex inside mesh
+    bool refineNeighboringFaces() const;
+
     void setZValueType( QgsMeshEditDigitizingAction::ZValueSource zValueSource );
 
   private slots:
     void updateSettings();
 
   private:
     QComboBox *mComboZValueType = nullptr;
+    QCheckBox *mCheckBoxRefineNeighboringFaces = nullptr;
+
+    friend class TestQgsMapToolEditMesh;
 };
 
 class APP_EXPORT QgsMapToolEditMeshFrame : public QgsMapToolAdvancedDigitizing

src/core/mesh/qgsmesheditor.cpp

@@ -640,6 +640,15 @@ QgsMeshEditingError QgsMeshEditor::removeFaces( const QList<int> &facesToRemove
   return error;
 }
 
+void QgsMeshEditor::addVertexWithDelaunayRefinement( const QgsMeshVertex &vertex, const double tolerance )
+{
+  int triangleIndex = mTriangularMesh->faceIndexForPoint_v2( vertex );
+  if ( triangleIndex == -1 )
+    return;
+
+  mUndoStack->push( new QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement( this, vertex, tolerance ) );
+}
+
 bool QgsMeshEditor::edgeCanBeFlipped( int vertexIndex1, int vertexIndex2 ) const
 {
   return mTopologicalMesh.edgeCanBeFlipped( vertexIndex1, vertexIndex2 );
@@ -1482,3 +1491,115 @@ void QgsMeshLayerUndoCommandAdvancedEditing::redo()
       mMeshEditor->applyEdit( edit );
   }
 }
+
+QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement::QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement(
+  QgsMeshEditor *meshEditor,
+  const QgsMeshVertex &vertex,
+  double tolerance )
+  : QgsMeshLayerUndoCommandMeshEdit( meshEditor )
+  , mVertex( vertex )
+  , mTolerance( tolerance )
+{
+  setText( QObject::tr( "Add vertex inside face with Delaunay refinement" ) );
+}
+
+void QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement::redo()
+{
+  if ( !mVertex.isEmpty() )
+  {
+    QgsMeshEditor::Edit edit;
+
+    mMeshEditor->applyAddVertex( edit, mVertex, mTolerance );
+    mEdits.append( edit );
+
+    QList<std::pair<int, int>> sharedEdges = innerEdges( secondNeighboringTriangularFaces() );
+
+    for ( std::pair<int, int> edge : sharedEdges )
+    {
+      if ( mMeshEditor->edgeCanBeFlipped( edge.first, edge.second ) && !mMeshEditor->topologicalMesh().delaunayConditionForEdge( edge.first, edge.second ) )
+      {
+        mMeshEditor->applyFlipEdge( edit, edge.first, edge.second );
+        mEdits.append( edit );
+      }
+    }
+
+    mVertex = QgsMeshVertex();
+  }
+  else
+  {
+    for ( QgsMeshEditor::Edit &edit : mEdits )
+      mMeshEditor->applyEdit( edit );
+  }
+}
+
+QSet<int> QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement::secondNeighboringTriangularFaces()
+{
+  const int vIndex = mMeshEditor->topologicalMesh().mesh()->vertexCount() - 1;
+  const QList<int> firstNeighborFaces = mMeshEditor->topologicalMesh().facesAroundVertex( vIndex );
+  QSet<int> firstNeighborVertices;
+  for ( int face : firstNeighborFaces )
+  {
+    const QgsMeshFace meshFace = mMeshEditor->topologicalMesh().mesh()->face( face );
+    for ( int vertex : meshFace )
+    {
+      firstNeighborVertices.insert( vertex );
+    }
+  }
+
+  QSet<int> secondNeighboringFaces;
+  for ( int vertex : firstNeighborVertices )
+  {
+    const QList<int> faces = mMeshEditor->topologicalMesh().facesAroundVertex( vertex );
+    for ( int face : faces )
+    {
+      if ( mMeshEditor->topologicalMesh().mesh()->face( face ).count() == 3 )
+        secondNeighboringFaces.insert( face );
+    }
+  }
+  return secondNeighboringFaces;
+}
+
+QList<std::pair<int, int>> QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement::innerEdges( const QSet<int> &faces )
+{
+  // edges and number of their occurrence in triangular faces
+  QMap<std::pair<int, int>, int> edges;
+
+  for ( int faceIndex : faces )
+  {
+    const QgsMeshFace face = mMeshEditor->topologicalMesh().mesh()->face( faceIndex );
+
+    for ( int i = 0; i < face.size(); i++ )
+    {
+      int next = i + 1;
+      if ( next == face.size() )
+      {
+        next = 0;
+      }
+
+      int minIndex = std::min( face.at( i ), face.at( next ) );
+      int maxIndex = std::max( face.at( i ), face.at( next ) );
+      std::pair<int, int> edge = std::pair<int, int>( minIndex, maxIndex );
+
+      int count = 1;
+      if ( edges.contains( edge ) )
+      {
+        count = edges.take( edge );
+        count++;
+      }
+
+      edges.insert( edge, count );
+    }
+  }
+
+  QList<std::pair<int, int>> sharedEdges;
+
+  for ( auto it = edges.begin(); it != edges.end(); it++ )
+  {
+    if ( it.value() == 2 )
+    {
+      sharedEdges.push_back( it.key() );
+    }
+  }
+
+  return sharedEdges;
+}

src/core/mesh/qgsmesheditor.h

@@ -305,6 +305,14 @@ class CORE_EXPORT QgsMeshEditor : public QObject
     //! Returns the maximum count of vertices per face that the mesh can support
     int maximumVerticesPerFace() const;
 
+    /**
+     * Add a vertex in a face with Delaunay refinement of neighboring faces
+     * All neighboring faces sharing a vertex will be refined to satisfy the Delaunay condition
+     *
+     * \since QGIS 3.42
+     */
+    void addVertexWithDelaunayRefinement( const QgsMeshVertex &vertex, const double tolerance );
+
   signals:
     //! Emitted when the mesh is edited
     void meshEdited();
@@ -363,6 +371,7 @@ class CORE_EXPORT QgsMeshEditor : public QObject
     friend class QgsMeshLayerUndoCommandFlipEdge;
     friend class QgsMeshLayerUndoCommandMerge;
     friend class QgsMeshLayerUndoCommandSplitFaces;
+    friend class QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement;
 
     friend class QgsMeshLayerUndoCommandAdvancedEditing;
 };
@@ -652,7 +661,28 @@ class QgsMeshLayerUndoCommandAdvancedEditing : public QgsMeshLayerUndoCommandMes
 };
 
 
+/**
+ * \ingroup core
+ *
+ * \brief  Class for undo/redo command for adding vertex to face with Delaunay Refiment of faces surrounding
+ *
+ * \since QGIS 3.42
+ */
+class QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement: public QgsMeshLayerUndoCommandMeshEdit
+{
+  public:
+
+    //! Constructor with the associated \a meshEditor and indexes \a vertex and \a tolerance
+    QgsMeshLayerUndoCommandAddVertexInFaceWithDelaunayRefinement( QgsMeshEditor *meshEditor, const QgsMeshVertex &vertex, double tolerance );
+
+    void redo() override;
+  private:
+    QList<std::pair<int, int>> innerEdges( const QSet<int> &faces );
+    QSet<int> secondNeighboringTriangularFaces();
 
+    QgsMeshVertex mVertex;
+    double mTolerance;
+};
 
 #endif //SIP_RUN
 

src/core/mesh/qgstopologicalmesh.cpp

@@ -19,6 +19,7 @@
 #include "qgsmesheditor.h"
 #include "qgsmessagelog.h"
 #include "qgsgeometryutils.h"
+#include "qgscircle.h"
 
 #include <poly2tri.h>
 #include <QSet>
@@ -2584,3 +2585,41 @@ QgsTopologicalMesh::Changes QgsTopologicalMesh::changeXYValue( const QList<int>
 
   return changes;
 }
+
+bool QgsTopologicalMesh::delaunayConditionForEdge( int vertexIndex1, int vertexIndex2 )
+{
+  int faceIndex1;
+  int faceIndex2;
+  int oppositeVertexFace1;
+  int oppositeVertexFace2;
+  int supposedOppositeVertexFace1;
+  int supposedoppositeVertexFace2;
+
+  bool result =  eitherSideFacesAndVertices(
+                   vertexIndex1,
+                   vertexIndex2,
+                   faceIndex1,
+                   faceIndex2,
+                   oppositeVertexFace1,
+                   supposedoppositeVertexFace2,
+                   supposedOppositeVertexFace1,
+                   oppositeVertexFace2 );
+
+  if ( ! result )
+    return false;
+
+  const QgsMeshFace face1 = mMesh->face( faceIndex1 );
+  const QgsMeshFace face2 = mMesh->face( faceIndex2 );
+
+  QgsCircle circle = QgsCircle::from3Points( mMesh->vertex( face1.at( 0 ) ),
+                     mMesh->vertex( face1.at( 1 ) ),
+                     mMesh->vertex( face1.at( 2 ) ) );
+  bool circle1ContainsPoint = circle.contains( mMesh->vertex( supposedoppositeVertexFace2 ) );
+
+  circle = QgsCircle::from3Points( mMesh->vertex( face2.at( 0 ) ),
+                                   mMesh->vertex( face2.at( 1 ) ),
+                                   mMesh->vertex( face2.at( 2 ) ) );
+  bool circle2ContainsPoint = circle.contains( mMesh->vertex( supposedOppositeVertexFace1 ) );
+
+  return !( circle1ContainsPoint || circle2ContainsPoint );
+}

src/core/mesh/qgstopologicalmesh.h

@@ -232,6 +232,14 @@ class CORE_EXPORT QgsTopologicalMesh
      */
     Changes flipEdge( int vertexIndex1, int vertexIndex2 );
 
+    /**
+     * Check if Delaunay condition holds for given edge
+     * returns TRUE if delaunay condition holds FALSE otherwise
+     *
+     * \since QGIS 3.42
+     */
+    bool delaunayConditionForEdge( int vertexIndex1, int vertexIndex2 );
+
     /**
      * Returns TRUE if faces separated by vertices with indexes \a vertexIndex1 and \a vertexIndex2 can be merged
      */