-
-
Notifications
You must be signed in to change notification settings - Fork 3k
/
DualEdgeTriangulation.h
249 lines (234 loc) · 14.6 KB
/
DualEdgeTriangulation.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
/***************************************************************************
DualEdgeTriangulation.h - description
-------------------
copyright : (C) 2004 by Marco Hugentobler
email : mhugent@geo.unizh.ch
***************************************************************************/
/***************************************************************************
* *
* 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. *
* *
***************************************************************************/
#ifndef DUALEDGETRIANGULATION_H
#define DUALEDGETRIANGULATION_H
#include "Triangulation.h"
#include "HalfEdge.h"
#include <QVector>
#include <QList>
#include <QSet>
#include <QColor>
#include <QFile>
#include <QTextStream>
#include <QMessageBox>
#include <cfloat>
#include <QBuffer>
#include <QStringList>
#include <QProgressDialog>
#include <QCursor>
#include "qgis_analysis.h"
/** \ingroup analysis
* DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data structure*/
class ANALYSIS_EXPORT DualEdgeTriangulation: public Triangulation
{
public:
DualEdgeTriangulation();
DualEdgeTriangulation( int nop, Triangulation *decorator );
virtual ~DualEdgeTriangulation();
void setDecorator( Triangulation *d ) {mDecorator = d;}
//! Adds a line (e.g. a break-, structure- or an isoline) to the triangulation. The class takes ownership of the line object and its points
void addLine( Line3D *line, bool breakline ) override;
//! Adds a point to the triangulation and returns the number of this point in case of success or -100 in case of failure
int addPoint( Point3D *p ) override;
//! Performs a consistency check, remove this later
virtual void performConsistencyTest() override;
//! Calculates the normal at a point on the surface
virtual bool calcNormal( double x, double y, Vector3D *result ) override;
//! Calculates x-, y and z-value of the point on the surface
virtual bool calcPoint( double x, double y, Point3D *result ) override;
//! Draws the points, edges and the forced lines
//virtual void draw(QPainter* p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height) const;
//! Returns a pointer to the point with number i
virtual Point3D *getPoint( unsigned int i ) const override;
//! Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedge)
int getOppositePoint( int p1, int p2 ) override;
//! Finds out, in which triangle the point with coordinates x and y is and assigns the numbers of the vertices to 'n1', 'n2' and 'n3' and the vertices to 'p1', 'p2' and 'p3'
//! \note not available in Python bindings
virtual bool getTriangle( double x, double y, Point3D *p1, int *n1, Point3D *p2, int *n2, Point3D *p3, int *n3 ) override;
//! Finds out, in which triangle the point with coordinates x and y is and assigns addresses to the points at the vertices to 'p1', 'p2' and 'p3
virtual bool getTriangle( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3 ) override;
//! Returns a pointer to a value list with the information of the triangles surrounding (counterclockwise) a point. Four integer values describe a triangle, the first three are the number of the half edges of the triangle and the fourth is -10, if the third (and most counterclockwise) edge is a breakline, and -20 otherwise. The value list has to be deleted by the code which called the method
QList<int> *getSurroundingTriangles( int pointno ) override;
//! Returns the largest x-coordinate value of the bounding box
virtual double getXMax() const override { return xMax; }
//! Returns the smallest x-coordinate value of the bounding box
virtual double getXMin() const override { return xMin; }
//! Returns the largest y-coordinate value of the bounding box
virtual double getYMax() const override { return yMax; }
//! Returns the smallest x-coordinate value of the bounding box
virtual double getYMin() const override { return yMin; }
//! Returns the number of points
virtual int getNumberOfPoints() const override;
//! Sets the behavior of the triangulation in case of crossing forced lines
virtual void setForcedCrossBehavior( Triangulation::ForcedCrossBehavior b ) override;
//! Sets the color of the normal edges
virtual void setEdgeColor( int r, int g, int b ) override;
//! Sets the color of the forced edges
virtual void setForcedEdgeColor( int r, int g, int b ) override;
//! Sets the color of the breaklines
virtual void setBreakEdgeColor( int r, int g, int b ) override;
//! Sets an interpolator object
void setTriangleInterpolator( TriangleInterpolator *interpolator ) override;
//! Eliminates the horizontal triangles by swapping or by insertion of new points
void eliminateHorizontalTriangles() override;
//! Adds points to make the triangles better shaped (algorithm of ruppert)
virtual void ruppertRefinement() override;
//! Returns true, if the point with coordinates x and y is inside the convex hull and false otherwise
bool pointInside( double x, double y ) override;
//! Reads the dual edge structure of a taff file
//bool readFromTAFF(QString fileName);
//! Saves the dual edge structure to a taff file
//bool saveToTAFF(QString fileName) const;
//! Swaps the edge which is closest to the point with x and y coordinates (if this is possible)
virtual bool swapEdge( double x, double y ) override;
//! 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
virtual QList<int> *getPointsAroundEdge( double x, double y ) override;
/** Saves the triangulation as a (line) shapefile
\return true in case of success*/
virtual bool saveAsShapefile( const QString &fileName ) const override;
protected:
//! X-coordinate of the upper right corner of the bounding box
double xMax;
//! X-coordinate of the lower left corner of the bounding box
double xMin;
//! Y-coordinate of the upper right corner of the bounding box
double yMax;
//! Y-coordinate of the lower left corner of the bounding box
double yMin;
//! Default value for the number of storable points at the beginning
static const unsigned int DEFAULT_STORAGE_FOR_POINTS = 100000;
//! Stores pointers to all points in the triangulations (including the points contained in the lines)
QVector<Point3D *> mPointVector;
//! Default value for the number of storable HalfEdges at the beginning
static const unsigned int DEFAULT_STORAGE_FOR_HALF_EDGES = 300006;
//! Stores pointers to the HalfEdges
QVector<HalfEdge *> mHalfEdge;
//! Association to an interpolator object
TriangleInterpolator *mTriangleInterpolator = nullptr;
//! Member to store the behavior in case of crossing forced segments
Triangulation::ForcedCrossBehavior mForcedCrossBehavior;
//! Color to paint the normal edges
QColor mEdgeColor;
//! Color to paint the forced edges
QColor mForcedEdgeColor;
//! Color to paint the breaklines
QColor mBreakEdgeColor;
//! Pointer to the decorator using this triangulation. It it is used directly, mDecorator equals this
Triangulation *mDecorator = nullptr;
//! Inserts an edge and makes sure, everything is ok with the storage of the edge. The number of the HalfEdge is returned
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, bool breakline );
//! 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
int baseEdgeOfTriangle( Point3D *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 );
//! Swaps 'edge' and test recursively for other swaps (delaunay criterion)
void doSwap( 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;
//! 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;
//! If an inserted point is exactly on an existing edge, 'baseEdgeOfTriangle' returns -20 and sets the variable 'mEdgeWithPoint'
unsigned int mEdgeWithPoint;
//! If an instability occurs in 'baseEdgeOfTriangle', mUnstableEdge is set to the value of the current edge
unsigned int mUnstableEdge;
//! If a point has been inserted twice, its number is stored in this member
int mTwiceInsPoint;
//! Returns true, if it is possible to swap an edge, otherwise false(concave quad or edge on (or outside) the convex hull)
bool swapPossible( unsigned int edge );
//! 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
void triangulatePolygon( QList<int> *poly, QList<int> *free, int mainedge );
//! 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
bool halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const;
//! Calculates the minimum angle, which would be present, if the specified halfedge would be swapped
double swapMinAngle( int edge ) const;
//! 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
int splitHalfEdge( int edge, float position );
//! Returns true, if a half edge is on the convex hull and false otherwise
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( Point3D *point, int edge, QSet<int> &set );
};
inline DualEdgeTriangulation::DualEdgeTriangulation()
: xMax( 0 )
, xMin( 0 )
, yMax( 0 )
, yMin( 0 )
, mTriangleInterpolator( nullptr )
, mForcedCrossBehavior( Triangulation::DeleteFirst )
, mEdgeColor( 0, 255, 0 )
, mForcedEdgeColor( 0, 0, 255 )
, mBreakEdgeColor( 100, 100, 0 )
, mDecorator( this )
, mEdgeInside( 0 )
, mEdgeOutside( 0 )
, mEdgeWithPoint( 0 )
, mUnstableEdge( 0 )
, mTwiceInsPoint( 0 )
{
mPointVector.reserve( DEFAULT_STORAGE_FOR_POINTS );
mHalfEdge.reserve( DEFAULT_STORAGE_FOR_HALF_EDGES );
}
inline DualEdgeTriangulation::DualEdgeTriangulation( int nop, Triangulation *decorator )
: xMax( 0 )
, xMin( 0 )
, yMax( 0 )
, yMin( 0 )
, mTriangleInterpolator( nullptr )
, mForcedCrossBehavior( Triangulation::DeleteFirst )
, mEdgeColor( 0, 255, 0 )
, mForcedEdgeColor( 0, 0, 255 )
, mBreakEdgeColor( 100, 100, 0 )
, mDecorator( decorator ? decorator : this )
, mEdgeInside( 0 )
, mEdgeOutside( 0 )
, mEdgeWithPoint( 0 )
, mUnstableEdge( 0 )
, mTwiceInsPoint( 0 )
{
mPointVector.reserve( nop );
mHalfEdge.reserve( nop );
}
inline int DualEdgeTriangulation::getNumberOfPoints() const
{
return mPointVector.count();
}
inline Point3D *DualEdgeTriangulation::getPoint( unsigned int i ) const
{
return mPointVector.at( i );
}
inline bool DualEdgeTriangulation::halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const
{
return (
( getPoint( mHalfEdge[edge]->getPoint() )->getX() >= xlowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->getX() <= xupright &&
getPoint( mHalfEdge[edge]->getPoint() )->getY() >= ylowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->getY() <= yupright ) ||
( getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getX() >= xlowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getX() <= xupright &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getY() >= ylowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getY() <= yupright )
);
}
#endif