fix QgsDuaEdgeTriangulation

.
qgis · Aug 9, 2020 · 32ebbd8 · 32ebbd8
1 parent 05b9482
commit 32ebbd8
Show file tree

Hide file tree

Showing 7 changed files with 613 additions and 254 deletions.
diff --git a/src/analysis/interpolation/qgsdualedgetriangulation.cpp b/src/analysis/interpolation/qgsdualedgetriangulation.cpp
diff --git a/src/analysis/interpolation/qgsdualedgetriangulation.h b/src/analysis/interpolation/qgsdualedgetriangulation.h
@@ -97,7 +97,7 @@ class ANALYSIS_EXPORT QgsDualEdgeTriangulation: public QgsTriangulation
 
     virtual QgsMesh triangulationToMesh() const override;
 
-  protected:
+  private:
     //! X-coordinate of the upper right corner of the bounding box
     double mXMax = 0;
     //! X-coordinate of the lower left corner of the bounding box
@@ -122,24 +122,31 @@ class ANALYSIS_EXPORT QgsDualEdgeTriangulation: public QgsTriangulation
     unsigned int insertEdge( int dual, int next, int point, bool mbreak, bool forced );
     //! Inserts a forced segment between the points with the numbers p1 and p2 into the triangulation and returns the number of a HalfEdge belonging to this forced edge or -100 in case of failure
     int insertForcedSegment( int p1, int p2, QgsInterpolator::SourceType segmentType );
-    //! Threshold for the leftOfTest to handle numerical instabilities
-    //const static double leftOfTresh=0.00001;
     //! Security to prevent endless loops in 'baseEdgeOfTriangle'. It there are more iteration then this number, the point will not be inserted
     static const int MAX_BASE_ITERATIONS = 300000;
     //! Returns the number of an edge which points to the point with number 'point' or -1 if there is an error
     int baseEdgeOfPoint( int point );
-    //! Returns the number of a HalfEdge from a triangle in which 'point' is in. If the number -10 is returned, this means, that 'point' is outside the convex hull. If -5 is returned, then numerical problems with the leftOfTest occurred (and the value of the possible edge is stored in the variable 'mUnstableEdge'. -20 means, that the inserted point is exactly on an edge (the number is stored in the variable 'mEdgeWithPoint'). -25 means, that the point is already in the triangulation (the number of the point is stored in the member 'mTwiceInsPoint'. If -100 is returned, this means that something else went wrong
+
+    /**
+     * Returns the number of a HalfEdge from a triangle in which 'point' is in. If the number -10 is returned, this means, that 'point' is outside the convex hull.
+     * If -5 is returned, then numerical problems with the leftOfTest occurred (and the value of the possible edge is stored in the variable 'mUnstableEdge'.
+     * -20 means, that the inserted point is exactly on an edge (the number is stored in the variable 'mEdgeWithPoint').
+     * -25 means, that the point is already in the triangulation (the number of the point is stored in the member 'mTwiceInsPoint'.
+     * If -100 is returned, this means that something else went wrong
+     */
     int baseEdgeOfTriangle( const QgsPoint &point );
     //! Checks, if 'edge' has to be swapped because of the empty circle criterion. If so, doSwap(...) is called.
-    bool checkSwap( unsigned int edge, unsigned int recursiveDeep );
+    bool checkSwapRecursivly( unsigned int edge, unsigned int recursiveDeep );
+    //! Return true if edge need to be swapped after Delaunay criteria control
+    bool isEdgeNeedSwap( unsigned int edge );
     //! Swaps 'edge' and test recursively for other swaps (delaunay criterion)
-    void doSwap( unsigned int edge, unsigned int recursiveDeep );
+    void doSwapRecursivly( unsigned int edge, unsigned int recursiveDeep );
     //! Swaps 'edge' and does no recursiv testing
     void doOnlySwap( unsigned int edge );
     //! Number of an edge which does not point to the virtual point. It continuously updated for a fast search
     unsigned int mEdgeInside = 0;
     //! Number of an edge on the outside of the convex hull. It is updated in method 'baseEdgeOfTriangle' to enable insertion of points outside the convex hull
-    unsigned int mEdgeOutside = 0;
+    int mEdgeOutside = -1;
     //! If an inserted point is exactly on an existing edge, 'baseEdgeOfTriangle' returns -20 and sets the variable 'mEdgeWithPoint'
     unsigned int mEdgeWithPoint = 0;
     //! If an instability occurs in 'baseEdgeOfTriangle', mUnstableEdge is set to the value of the current edge
@@ -160,6 +167,11 @@ class ANALYSIS_EXPORT QgsDualEdgeTriangulation: public QgsTriangulation
     bool edgeOnConvexHull( int edge );
     //! Function needed for the ruppert algorithm. Tests, if point is in the circle through both endpoints of edge and the endpoint of edge->dual->next->point. If so, the function calls itself recursively for edge->next and edge->next->next. Stops, if it finds a forced edge or a convex hull edge
     void evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set );
+    //! Dimension of the triangulation, -1 : no point, 0 : 1 point, 1 : 2 point or only colinear edges, 2 : 3 or more points and at least 2 non colinear edges
+    int mDimension = -1;
+
+    int firstEdgeOutSide();
+
 
     friend class TestQgsInterpolator;
 };

diff --git a/src/core/mesh/qgsmeshdataprovider.cpp b/src/core/mesh/qgsmeshdataprovider.cpp
@@ -131,6 +131,26 @@ void QgsMesh::clear()
   faces.clear();
 }
 
+bool QgsMesh::compareFaces( const QgsMeshFace &face1, const QgsMeshFace &face2 )
+{
+  if ( face1.count() != face2.count() )
+    return false;
+
+  int startFace2 = 0;
+  for ( int i = 0; i < face2.count(); ++i )
+    if ( face2.at( i ) == face1.at( 0 ) )
+    {
+      startFace2 = i;
+      break;
+    }
+
+  for ( int i = 0; i < face1.count(); ++i )
+    if ( face1.at( i ) != face2.at( ( i + startFace2 ) % ( face2.count() ) ) )
+      return false;
+
+  return true;
+}
+
 bool QgsMesh::contains( const QgsMesh::ElementType &type ) const
 {
   switch ( type )

diff --git a/src/core/mesh/qgsmeshdataprovider.h b/src/core/mesh/qgsmeshdataprovider.h
@@ -102,6 +102,12 @@ struct CORE_EXPORT QgsMesh
     */
   void clear();
 
+  /**
+   * Compare two faces, return true if they are equivalent : same indexes and same clock wise
+    * \since QGIS 3.16
+   */
+  static bool compareFaces( const QgsMeshFace &face1, const QgsMeshFace &face2 );
+
   QVector<QgsMeshVertex> vertices SIP_SKIP;
   QVector<QgsMeshEdge> edges SIP_SKIP;
   QVector<QgsMeshFace> faces SIP_SKIP;

diff --git a/tests/src/analysis/CMakeLists.txt b/tests/src/analysis/CMakeLists.txt
@@ -92,6 +92,7 @@ SET(TESTS
  testqgsninecellfilters.cpp 
  testqgsmeshcalculator.cpp
  testqgsmeshcontours.cpp
+ testqgstriangulation.cpp
 )
 
 FOREACH(TESTSRC ${TESTS})

diff --git a/tests/src/analysis/testqgsinterpolator.cpp b/tests/src/analysis/testqgsinterpolator.cpp
@@ -274,11 +274,12 @@ void TestQgsInterpolator::dualEdge()
   QCOMPARE( mesh.vertex( 3 ), QgsMeshVertex( 2.0, 4.0, 6.0 ) );
   QCOMPARE( mesh.vertex( 4 ), QgsMeshVertex( 2.0, 2.0, 7.0 ) );
 
-  QCOMPARE( mesh.face( 0 ), QgsMeshFace( {0, 4, 3} ) );
-  QCOMPARE( mesh.face( 1 ), QgsMeshFace( {1, 4, 0} ) );
-  QCOMPARE( mesh.face( 2 ), QgsMeshFace( {2, 4, 1} ) );
-  QCOMPARE( mesh.face( 3 ), QgsMeshFace( {4, 2, 3} ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 4, 3} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 1 ), QgsMeshFace( {4, 2, 3} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 2 ), QgsMeshFace( {1, 4, 0} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 3 ), QgsMeshFace( {2, 4, 1} ) ) );
 }
 
+
 QGSTEST_MAIN( TestQgsInterpolator )
 #include "testqgsinterpolator.moc"
diff --git a/tests/src/analysis/testqgstriangulation.cpp b/tests/src/analysis/testqgstriangulation.cpp
@@ -0,0 +1,156 @@
+/***************************************************************************
+  testqgsinterpolator.cpp
+  -----------------------
+Date                 : August 2020
+Copyright            : (C) 2020 by Vincent Cloarec
+Email                : vcloarec 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#include "qgstest.h"
+
+#include "qgsapplication.h"
+#include "qgsdualedgetriangulation.h"
+
+class TestQgsTriangulation : public QObject
+{
+    Q_OBJECT
+
+  public:
+
+  private slots:
+    void initTestCase();// will be called before the first testfunction is executed.
+    void cleanupTestCase();// will be called after the last testfunction was executed.
+    void init() ;// will be called before each testfunction is executed.
+    void cleanup() ;// will be called after every testfunction.
+    void dualEdge();
+
+  private:
+};
+
+void  TestQgsTriangulation::initTestCase()
+{}
+
+void TestQgsTriangulation::cleanupTestCase()
+{}
+
+void TestQgsTriangulation::init()
+{}
+
+void TestQgsTriangulation::cleanup()
+{}
+
+void TestQgsTriangulation::dualEdge()
+{
+  //3 points
+  QgsDualEdgeTriangulation triangulation;
+  // Add colinear points
+  triangulation.addPoint( QgsPoint( 1, 0, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 1, 0 ) );
+  QgsMesh mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 2 );
+  QCOMPARE( mesh.faceCount(), 0 );
+  triangulation.addPoint( QgsPoint( 2, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 3 );
+  QCOMPARE( mesh.faceCount(), 1 );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 2, 1} ) ) );
+
+  //4 points
+  triangulation = QgsDualEdgeTriangulation();
+  // Add colinear points
+  triangulation.addPoint( QgsPoint( 1, 0, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 1, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 2 );
+  QCOMPARE( mesh.faceCount(), 0 );
+  triangulation.addPoint( QgsPoint( 2, 2, 0 ) );
+  triangulation.addPoint( QgsPoint( 2, 3, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 4 );
+  QCOMPARE( mesh.faceCount(), 2 );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 2, 1} ) ) );
+
+  //3 first colinear points
+  triangulation = QgsDualEdgeTriangulation();
+  triangulation.addPoint( QgsPoint( 1, 0, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 1, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 3 );
+  QCOMPARE( mesh.faceCount(), 0 );
+  triangulation.addPoint( QgsPoint( 2, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 4 );
+  QCOMPARE( mesh.faceCount(), 2 );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 3, 1} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 1 ), QgsMeshFace( {3, 2, 1} ) ) );
+  triangulation.addPoint( QgsPoint( 2, 3, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 5 );
+  QCOMPARE( mesh.faceCount(), 3 );
+
+  //3 first colinear points with different order
+  triangulation = QgsDualEdgeTriangulation();
+  triangulation.addPoint( QgsPoint( 1, 0, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 2, 0 ) );
+  triangulation.addPoint( QgsPoint( 1, 1, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 3 );
+  QCOMPARE( mesh.faceCount(), 0 );
+  triangulation.addPoint( QgsPoint( 2, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 4 );
+  QCOMPARE( mesh.faceCount(), 2 );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 3, 2} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 1 ), QgsMeshFace( {1, 2, 3} ) ) );
+  triangulation.addPoint( QgsPoint( 2, 3, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 5 );
+  QCOMPARE( mesh.faceCount(), 3 );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 0 ), QgsMeshFace( {0, 3, 2} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 1 ), QgsMeshFace( {1, 3, 4} ) ) );
+  QVERIFY( QgsMesh::compareFaces( mesh.face( 2 ), QgsMeshFace( {1, 2, 3} ) ) );
+
+  //4 first colinear points
+  triangulation = QgsDualEdgeTriangulation();
+  triangulation.addPoint( QgsPoint( 1, 1, 0 ) );
+  triangulation.addPoint( QgsPoint( 2, 1, 0 ) );
+  triangulation.addPoint( QgsPoint( 3, 1, 0 ) );
+  triangulation.addPoint( QgsPoint( 4, 1, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 4 );
+  QCOMPARE( mesh.faceCount(), 0 );
+  triangulation.addPoint( QgsPoint( 1, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 5 );
+  QCOMPARE( mesh.faceCount(), 3 );
+  triangulation.addPoint( QgsPoint( 2, 2, 0 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.vertices.count(), 6 );
+  QCOMPARE( mesh.faceCount(), 4 );
+
+  triangulation = QgsDualEdgeTriangulation();
+  triangulation.addPoint( QgsPoint( 2, 0, 1 ) );
+  triangulation.addPoint( QgsPoint( 0, 2, 1 ) );
+  triangulation.addPoint( QgsPoint( 2, 4, 1 ) );
+  triangulation.addPoint( QgsPoint( 4, 2, 1 ) );
+
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.faceCount(), 2 );
+  QCOMPARE( mesh.vertexCount(), 4 );
+
+  //add point exactly on existing edge
+  triangulation.addPoint( QgsPoint( 2, 2, 1 ) );
+  mesh = triangulation.triangulationToMesh();
+  QCOMPARE( mesh.faceCount(), 4 );
+  QCOMPARE( mesh.vertexCount(), 5 );
+}
+
+QGSTEST_MAIN( TestQgsTriangulation )
+#include "testqgstriangulation.moc"