Branch data Line data Source code
1 : : /**************************************************************************
2 : : QgsDualEdgeTriangulation.h - description
3 : : -------------------
4 : : copyright : (C) 2004 by Marco Hugentobler
5 : : email : mhugent@geo.unizh.ch
6 : : ***************************************************************************/
7 : :
8 : : /***************************************************************************
9 : : * *
10 : : * This program is free software; you can redistribute it and/or modify *
11 : : * it under the terms of the GNU General Public License as published by *
12 : : * the Free Software Foundation; either version 2 of the License, or *
13 : : * (at your option) any later version. *
14 : : * *
15 : : ***************************************************************************/
16 : :
17 : : #ifndef DUALEDGETRIANGULATION_H
18 : : #define DUALEDGETRIANGULATION_H
19 : :
20 : : #include <QVector>
21 : : #include <QList>
22 : : #include <QSet>
23 : : #include <QColor>
24 : : #include <QFile>
25 : : #include <QTextStream>
26 : : #include <QMessageBox>
27 : : #include <QBuffer>
28 : : #include <QStringList>
29 : : #include <QCursor>
30 : :
31 : : #include <cfloat>
32 : :
33 : : #include "qgis_sip.h"
34 : : #include "qgis_analysis.h"
35 : : #include "qgspoint.h"
36 : :
37 : : #include "qgstriangulation.h"
38 : : #include "HalfEdge.h"
39 : :
40 : : #define SIP_NO_FILE
41 : :
42 : : /**
43 : : * \ingroup analysis
44 : : * \brief DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data structure.
45 : : * \note Not available in Python bindings.
46 : : *
47 : : * \since QGIS 3.16
48 : : */
49 : : class ANALYSIS_EXPORT QgsDualEdgeTriangulation: public QgsTriangulation
50 : : {
51 : : public:
52 : : QgsDualEdgeTriangulation();
53 : : //! Constructor with a number of points to reserve
54 : : QgsDualEdgeTriangulation( int nop );
55 : : ~QgsDualEdgeTriangulation() override;
56 : : void addLine( const QVector< QgsPoint > &points, QgsInterpolator::SourceType lineType ) override;
57 : : int addPoint( const QgsPoint &p ) override;
58 : : //! Performs a consistency check, remove this later
59 : : void performConsistencyTest() override;
60 : : //! Calculates the normal at a point on the surface
61 : : bool calcNormal( double x, double y, QgsPoint &result SIP_OUT ) override;
62 : : bool calcPoint( double x, double y, QgsPoint &result SIP_OUT ) override;
63 : : //! Draws the points, edges and the forced lines
64 : : //virtual void draw(QPainter* p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height) const;
65 : : //! Returns a pointer to the point with number i
66 : : QgsPoint *point( int i ) const override;
67 : : int oppositePoint( int p1, int p2 ) override;
68 : : bool triangleVertices( double x, double y, QgsPoint &p1 SIP_OUT, int &n1 SIP_OUT, QgsPoint &p2 SIP_OUT, int &n2 SIP_OUT, QgsPoint &p3 SIP_OUT, int &n3 SIP_OUT ) override;
69 : : bool triangleVertices( double x, double y, QgsPoint &p1 SIP_OUT, QgsPoint &p2 SIP_OUT, QgsPoint &p3 SIP_OUT ) override;
70 : : QList<int> surroundingTriangles( int pointno ) override;
71 : : //! Returns the largest x-coordinate value of the bounding box
72 : 0 : double xMax() const override { return mXMax; }
73 : : //! Returns the smallest x-coordinate value of the bounding box
74 : 0 : double xMin() const override { return mXMin; }
75 : : //! Returns the largest y-coordinate value of the bounding box
76 : 0 : double yMax() const override { return mYMax; }
77 : : //! Returns the smallest x-coordinate value of the bounding box
78 : 0 : double yMin() const override { return mYMin; }
79 : : //! Returns the number of points
80 : : int pointsCount() const override;
81 : : //! Sets the behavior of the triangulation in case of crossing forced lines
82 : : void setForcedCrossBehavior( QgsTriangulation::ForcedCrossBehavior b ) override;
83 : : //! Sets an interpolator object
84 : : void setTriangleInterpolator( TriangleInterpolator *interpolator ) override;
85 : : //! Eliminates the horizontal triangles by swapping or by insertion of new points
86 : : void eliminateHorizontalTriangles() override;
87 : : //! Adds points to make the triangles better shaped (algorithm of ruppert)
88 : : void ruppertRefinement() override;
89 : : //! Returns TRUE, if the point with coordinates x and y is inside the convex hull and FALSE otherwise
90 : : bool pointInside( double x, double y ) override;
91 : : //! Swaps the edge which is closest to the point with x and y coordinates (if this is possible)
92 : : bool swapEdge( double x, double y ) override;
93 : : //! Returns a value list with the numbers of the four points, which would be affected by an edge swap. This function is e.g. needed by NormVecDecorator to know the points, for which the normals have to be recalculated. The returned ValueList has to be deleted by the code which calls the method
94 : : QList<int> pointsAroundEdge( double x, double y ) override;
95 : :
96 : : bool saveTriangulation( QgsFeatureSink *sink, QgsFeedback *feedback = nullptr ) const override;
97 : :
98 : : virtual QgsMesh triangulationToMesh( QgsFeedback *feedback = nullptr ) const override;
99 : :
100 : : private:
101 : : //! X-coordinate of the upper right corner of the bounding box
102 : 0 : double mXMax = 0;
103 : : //! X-coordinate of the lower left corner of the bounding box
104 : 0 : double mXMin = 0;
105 : : //! Y-coordinate of the upper right corner of the bounding box
106 : 0 : double mYMax = 0;
107 : : //! Y-coordinate of the lower left corner of the bounding box
108 : 0 : double mYMin = 0;
109 : : //! Default value for the number of storable points at the beginning
110 : : static const unsigned int DEFAULT_STORAGE_FOR_POINTS = 100000;
111 : : //! Stores pointers to all points in the triangulations (including the points contained in the lines)
112 : : QVector<QgsPoint *> mPointVector;
113 : : //! Default value for the number of storable HalfEdges at the beginning
114 : : static const unsigned int DEFAULT_STORAGE_FOR_HALF_EDGES = 300006;
115 : : //! Stores pointers to the HalfEdges
116 : : QVector<HalfEdge *> mHalfEdge;
117 : : //! Association to an interpolator object
118 : 0 : TriangleInterpolator *mTriangleInterpolator = nullptr;
119 : : //! Member to store the behavior in case of crossing forced segments
120 : 0 : QgsTriangulation::ForcedCrossBehavior mForcedCrossBehavior = QgsTriangulation::DeleteFirst;
121 : : //! Inserts an edge and makes sure, everything is OK with the storage of the edge. The number of the HalfEdge is returned
122 : : unsigned int insertEdge( int dual, int next, int point, bool mbreak, bool forced );
123 : : //! 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
124 : : int insertForcedSegment( int p1, int p2, QgsInterpolator::SourceType segmentType );
125 : : //! Security to prevent endless loops in 'baseEdgeOfTriangle'. It there are more iteration then this number, the point will not be inserted
126 : : static const int MAX_BASE_ITERATIONS = 300000;
127 : : //! Returns the number of an edge which points to the point with number 'point' or -1 if there is an error
128 : : int baseEdgeOfPoint( int point );
129 : :
130 : : /**
131 : : * 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.
132 : : * If -5 is returned, then numerical problems with the leftOfTest occurred (and the value of the possible edge is stored in the variable 'mUnstableEdge'.
133 : : * -20 means, that the inserted point is exactly on an edge (the number is stored in the variable 'mEdgeWithPoint').
134 : : * -25 means, that the point is already in the triangulation (the number of the point is stored in the member 'mTwiceInsPoint'.
135 : : * If -100 is returned, this means that something else went wrong
136 : : */
137 : : int baseEdgeOfTriangle( const QgsPoint &point );
138 : : //! Checks, if 'edge' has to be swapped because of the empty circle criterion. If so, doSwap(...) is called.
139 : : bool checkSwapRecursively( unsigned int edge, unsigned int recursiveDeep );
140 : : //! Return true if edge need to be swapped after Delaunay criteria control
141 : : bool isEdgeNeedSwap( unsigned int edge ) const;
142 : : //! Swaps 'edge' and test recursively for other swaps (delaunay criterion)
143 : : void doSwapRecursively( unsigned int edge, unsigned int recursiveDeep );
144 : : //! Swaps 'edge' and does no recursiv testing
145 : : void doOnlySwap( unsigned int edge );
146 : : //! Number of an edge which does not point to the virtual point. It continuously updated for a fast search
147 : 0 : unsigned int mEdgeInside = 0;
148 : : //! 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
149 : 0 : int mEdgeOutside = -1;
150 : : //! If an inserted point is exactly on an existing edge, 'baseEdgeOfTriangle' returns -20 and sets the variable 'mEdgeWithPoint'
151 : 0 : unsigned int mEdgeWithPoint = 0;
152 : : //! If an instability occurs in 'baseEdgeOfTriangle', mUnstableEdge is set to the value of the current edge
153 : 0 : unsigned int mUnstableEdge = 0;
154 : : //! If a point has been inserted twice, its number is stored in this member
155 : 0 : int mTwiceInsPoint = 0;
156 : : //! Returns TRUE, if it is possible to swap an edge, otherwise FALSE(concave quad or edge on (or outside) the convex hull)
157 : : bool swapPossible( unsigned int edge ) const;
158 : : //! Divides a polygon in a triangle and two polygons and calls itself recursively for these two polygons. 'poly' is a pointer to a list with the numbers of the edges of the polygon, 'free' is a pointer to a list of free halfedges, and 'mainedge' is the number of the edge, towards which the new triangle is inserted. Mainedge has to be the same as poly->begin(), otherwise the recursion does not work
159 : : void triangulatePolygon( QList<int> *poly, QList<int> *free, int mainedge );
160 : : //! Tests, if the bounding box of the halfedge with index i intersects the specified bounding box. The main purpose for this method is the drawing of the triangulation
161 : : bool halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const;
162 : : //! Calculates the minimum angle, which would be present, if the specified halfedge would be swapped
163 : : double swapMinAngle( int edge ) const;
164 : : //! Inserts a new point on the halfedge with number 'edge'. The position can have a value from 0 to 1 (e.g. 0.5 would be in the middle). The return value is the number of the new inserted point. tin is the triangulation, which should be used to calculate the elevation of the inserted point
165 : : int splitHalfEdge( int edge, float position );
166 : : //! Returns TRUE, if a half edge is on the convex hull and FALSE otherwise
167 : : bool edgeOnConvexHull( int edge );
168 : : //! 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
169 : : void evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set );
170 : : //! 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
171 : 0 : int mDimension = -1;
172 : :
173 : : int firstEdgeOutSide();
174 : :
175 : : void removeLastPoint();
176 : :
177 : :
178 : : friend class TestQgsInterpolator;
179 : : };
180 : :
181 : : #ifndef SIP_RUN
182 : :
183 : 0 : inline QgsDualEdgeTriangulation::QgsDualEdgeTriangulation()
184 : 0 : {
185 : 0 : mPointVector.reserve( DEFAULT_STORAGE_FOR_POINTS );
186 : 0 : mHalfEdge.reserve( DEFAULT_STORAGE_FOR_HALF_EDGES );
187 : 0 : }
188 : :
189 : 0 : inline QgsDualEdgeTriangulation::QgsDualEdgeTriangulation( int nop )
190 : 0 : {
191 : 0 : mPointVector.reserve( nop );
192 : 0 : mHalfEdge.reserve( nop );
193 : 0 : }
194 : :
195 : 0 : inline int QgsDualEdgeTriangulation::pointsCount() const
196 : : {
197 : 0 : return mPointVector.count();
198 : : }
199 : :
200 : 0 : inline QgsPoint *QgsDualEdgeTriangulation::point( int i ) const
201 : : {
202 : 0 : if ( i < 0 || i >= mPointVector.count() )
203 : 0 : return nullptr;
204 : :
205 : 0 : return mPointVector.at( i );
206 : 0 : }
207 : :
208 : : inline bool QgsDualEdgeTriangulation::halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const
209 : : {
210 : : return (
211 : : ( point( mHalfEdge[edge]->getPoint() )->x() >= xlowleft &&
212 : : point( mHalfEdge[edge]->getPoint() )->x() <= xupright &&
213 : : point( mHalfEdge[edge]->getPoint() )->y() >= ylowleft &&
214 : : point( mHalfEdge[edge]->getPoint() )->y() <= yupright ) ||
215 : : ( point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->x() >= xlowleft &&
216 : : point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->x() <= xupright &&
217 : : point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->y() >= ylowleft &&
218 : : point( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->y() <= yupright )
219 : : );
220 : : }
221 : :
222 : : #endif
223 : : #endif
224 : :
225 : :
226 : :
|
