Skip to content

Commit 5d9b66f

Browse files
author
mhugent
committedJul 30, 2009
FEATURE: possibility to add line layers as constrains for triangulation in interpolation plugin
git-svn-id: http://svn.osgeo.org/qgis/trunk@11211 c8812cc2-4d05-0410-92ff-de0c093fc19c
1 parent 257e30e commit 5d9b66f

16 files changed

+500
-182
lines changed
 

‎src/plugins/interpolation/DualEdgeTriangulation.cc

Lines changed: 92 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,9 @@
1717

1818
#include "DualEdgeTriangulation.h"
1919
#include <map>
20+
#include "qgsgeometry.h"
2021
#include "qgslogger.h"
22+
#include "qgsvectorfilewriter.h"
2123

2224
double leftOfTresh = 0.00000001;
2325

@@ -75,7 +77,7 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline )
7577
for ( i = 0;i < line->getSize();i++ )
7678
{
7779
line->goToNext();
78-
actpoint = mDecorator->addPoint( line->getPoint() );
80+
actpoint = /*mDecorator->*/addPoint( line->getPoint() );
7981
if ( actpoint != -100 )
8082
{
8183
i++;
@@ -85,20 +87,22 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline )
8587

8688
if ( actpoint == -100 )//no point of the line could be inserted
8789
{
90+
delete line;
8891
return;
8992
}
9093

9194
for ( ;i < line->getSize();i++ )
9295
{
9396
line->goToNext();
94-
currentpoint = mDecorator->addPoint( line->getPoint() );
95-
//if(currentpoint!=-100&&actpoint!=-100&&currentpoint!=actpoint)//-100 is the return value if the point could not be not inserted
96-
//{
97-
// insertForcedSegment(actpoint,currentpoint,breakline);
98-
//}
97+
currentpoint = /*mDecorator->*/addPoint( line->getPoint() );
98+
if ( currentpoint != -100 && actpoint != -100 && currentpoint != actpoint )//-100 is the return value if the point could not be not inserted
99+
{
100+
insertForcedSegment( actpoint, currentpoint, breakline );
101+
}
99102
actpoint = currentpoint;
100103
}
101104
}
105+
delete line;
102106
}
103107

104108
int DualEdgeTriangulation::addPoint( Point3D* p )
@@ -1178,7 +1182,6 @@ unsigned int DualEdgeTriangulation::insertEdge( int dual, int next, int point, b
11781182

11791183
}
11801184

1181-
#if 0
11821185
int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
11831186
{
11841187
if ( p1 == p2 )
@@ -1436,8 +1439,8 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
14361439
}
14371440

14381441
//set the flags 'forced' and 'break' to false for every edge and dualedge of 'crossEdges'
1439-
QListIterator<int> iter;
1440-
for ( iter = crossedEdges.begin();iter != crossedEdges.end();++iter )
1442+
QList<int>::const_iterator iter;
1443+
for ( iter = crossedEdges.constBegin();iter != crossedEdges.constEnd();++iter )
14411444
{
14421445
mHalfEdge[( *( iter ) )]->setForced( false );
14431446
mHalfEdge[( *( iter ) )]->setBreak( false );
@@ -1467,8 +1470,9 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
14671470

14681471
//finish the polygon on the left side
14691472
int actpointl = p2;
1470-
QListIterator<int> leftiter;
1471-
leftiter = crossedEdges.fromLast();
1473+
QList<int>::const_iterator leftiter; //todo: is there a better way to set an iterator to the last list element?
1474+
leftiter = crossedEdges.constEnd();
1475+
--leftiter;
14721476
while ( true )
14731477
{
14741478
int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext()]->getPoint();
@@ -1479,7 +1483,7 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
14791483
int theedge = mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext();
14801484
leftPolygon.append( theedge );
14811485
}
1482-
if ( leftiter == crossedEdges.begin() )
1486+
if ( leftiter == crossedEdges.constBegin() )
14831487
{break;}
14841488
--leftiter;
14851489
}
@@ -1488,9 +1492,9 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
14881492
leftPolygon.append( mHalfEdge[crossedEdges.first()]->getNext() );
14891493

14901494
//finish the polygon on the right side
1491-
QValueListIterator<int> rightiter;
1495+
QList<int>::const_iterator rightiter;
14921496
int actpointr = p1;
1493-
for ( rightiter = crossedEdges.begin();rightiter != crossedEdges.end();++rightiter )
1497+
for ( rightiter = crossedEdges.constBegin();rightiter != crossedEdges.constEnd();++rightiter )
14941498
{
14951499
int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext()]->getPoint();
14961500
if ( newpoint != actpointr )
@@ -1509,15 +1513,17 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
15091513

15101514
//set the necessary nexts of leftPolygon(exept the first)
15111515
int actedgel = leftPolygon[1];
1512-
for ( leftiter = leftPolygon.at( 2 );leftiter != leftPolygon.end();++leftiter )
1516+
leftiter = leftPolygon.constBegin(); leftiter += 2;
1517+
for ( ;leftiter != leftPolygon.constEnd();++leftiter )
15131518
{
15141519
mHalfEdge[actedgel]->setNext(( *leftiter ) );
15151520
actedgel = ( *leftiter );
15161521
}
15171522

15181523
//set all the necessary nexts of rightPolygon
15191524
int actedger = rightPolygon[1];
1520-
for ( rightiter = rightPolygon.at( 2 );rightiter != rightPolygon.end();++rightiter )
1525+
rightiter = rightPolygon.constBegin(); rightiter += 2;
1526+
for ( ;rightiter != rightPolygon.constEnd();++rightiter )
15211527
{
15221528
mHalfEdge[actedger]->setNext(( *rightiter ) );
15231529
actedger = ( *( rightiter ) );
@@ -1532,8 +1538,6 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
15321538
mHalfEdge[rightPolygon.first()]->setPoint( p1 );
15331539
mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
15341540

1535-
1536-
15371541
triangulatePolygon( &leftPolygon, &freelist, firstedge );
15381542
triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );
15391543

@@ -1545,7 +1549,6 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
15451549

15461550
return leftPolygon.first();
15471551
}
1548-
#endif //0
15491552

15501553
void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b )
15511554
{
@@ -2442,7 +2445,6 @@ bool DualEdgeTriangulation::swapPossible( unsigned int edge )
24422445
return true;
24432446
}
24442447

2445-
#if 0
24462448
void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* free, int mainedge )
24472449
{
24482450
if ( poly && free )
@@ -2453,13 +2455,13 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
24532455
}
24542456

24552457
//search for the edge pointing on the closest point(distedge) and for the next(nextdistedge)
2456-
QValueListIterator<int> iterator = ++( poly->begin() );//go to the second edge
2458+
QList<int>::const_iterator iterator = ++( poly->constBegin() );//go to the second edge
24572459
double distance = MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
24582460
int distedge = ( *iterator );
24592461
int nextdistedge = mHalfEdge[( *iterator )]->getNext();
24602462
++iterator;
24612463

2462-
while ( iterator != --( poly->end() ) )
2464+
while ( iterator != --( poly->constEnd() ) )
24632465
{
24642466
if ( MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] ) < distance )
24652467
{
@@ -2476,15 +2478,15 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
24762478
int insertb = mHalfEdge[inserta]->getDual();
24772479
free->pop_front();
24782480

2479-
mHalfEdge[inserta]->setNext(( *( poly->at( 1 ) ) ) );
2481+
mHalfEdge[inserta]->setNext(( poly->at( 1 ) ) );
24802482
mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
24812483
mHalfEdge[insertb]->setNext( nextdistedge );
24822484
mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
24832485
mHalfEdge[distedge]->setNext( inserta );
24842486
mHalfEdge[mainedge]->setNext( insertb );
24852487

2486-
QValueList<int> polya;
2487-
for ( iterator = ( ++( poly->begin() ) );( *iterator ) != nextdistedge;++iterator )
2488+
QList<int> polya;
2489+
for ( iterator = ( ++( poly->constBegin() ) );( *iterator ) != nextdistedge;++iterator )
24882490
{
24892491
polya.append(( *iterator ) );
24902492
}
@@ -2507,16 +2509,16 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
25072509
int insertb = mHalfEdge[inserta]->getDual();
25082510
free->pop_front();
25092511

2510-
mHalfEdge[inserta]->setNext(( *( poly->at( 2 ) ) ) );
2512+
mHalfEdge[inserta]->setNext(( poly->at( 2 ) ) );
25112513
mHalfEdge[inserta]->setPoint( mHalfEdge[distedge]->getPoint() );
25122514
mHalfEdge[insertb]->setNext( mainedge );
25132515
mHalfEdge[insertb]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
25142516
mHalfEdge[distedge]->setNext( insertb );
25152517
mHalfEdge[( *( --poly->end() ) )]->setNext( inserta );
25162518

2517-
QValueList<int> polya;
2518-
iterator = poly->at( 2 );
2519-
while ( iterator != poly->end() )
2519+
QList<int> polya;
2520+
iterator = poly->constBegin(); iterator += 2;
2521+
while ( iterator != poly->constEnd() )
25202522
{
25212523
polya.append(( *iterator ) );
25222524
++iterator;
@@ -2536,7 +2538,7 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
25362538
int insertd = mHalfEdge[insertc]->getDual();
25372539
free->pop_front();
25382540

2539-
mHalfEdge[inserta]->setNext(( *( poly->at( 1 ) ) ) );
2541+
mHalfEdge[inserta]->setNext(( poly->at( 1 ) ) );
25402542
mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
25412543
mHalfEdge[insertb]->setNext( insertd );
25422544
mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
@@ -2550,17 +2552,17 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
25502552
mHalfEdge[( *( --poly->end() ) )]->setNext( insertc );
25512553

25522554
//build two new polygons for recursive triangulation
2553-
QValueList<int> polya;
2554-
QValueList<int> polyb;
2555+
QList<int> polya;
2556+
QList<int> polyb;
25552557

2556-
for ( iterator = ++( poly->begin() );( *iterator ) != nextdistedge;++iterator )
2558+
for ( iterator = ++( poly->constBegin() );( *iterator ) != nextdistedge;++iterator )
25572559
{
25582560
polya.append(( *iterator ) );
25592561
}
25602562
polya.prepend( inserta );
25612563

25622564

2563-
while ( iterator != poly->end() )
2565+
while ( iterator != poly->constEnd() )
25642566
{
25652567
polyb.append(( *iterator ) );
25662568
++iterator;
@@ -2578,7 +2580,6 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
25782580
}
25792581

25802582
}
2581-
#endif //0
25822583

25832584
bool DualEdgeTriangulation::pointInside( double x, double y )
25842585
{
@@ -3070,6 +3071,62 @@ QList<int>* DualEdgeTriangulation::getPointsAroundEdge( double x, double y )
30703071
}
30713072
}
30723073

3074+
bool DualEdgeTriangulation::saveAsShapefile( const QString& fileName ) const
3075+
{
3076+
QgsFieldMap fields;
3077+
fields.insert( 0, QgsField( "type", QVariant::String, "String" ) );
3078+
QgsVectorFileWriter writer( fileName, "Utf-8", fields, QGis::WKBLineString, 0 );
3079+
if ( writer.hasError() != QgsVectorFileWriter::NoError )
3080+
{
3081+
return false;
3082+
}
3083+
3084+
bool alreadyVisitedEdges[mHalfEdge.size()];
3085+
3086+
for ( int i = 0; i < mHalfEdge.size(); ++i )
3087+
{
3088+
alreadyVisitedEdges[i] = false;
3089+
}
3090+
3091+
for ( int i = 0; i < mHalfEdge.size(); ++i )
3092+
{
3093+
HalfEdge* currentEdge = mHalfEdge[i];
3094+
if ( currentEdge->getPoint() != -1 && mHalfEdge[currentEdge->getDual()]->getPoint() != -1 && !alreadyVisitedEdges[currentEdge->getDual()] )
3095+
{
3096+
QgsFeature edgeLineFeature;
3097+
3098+
//geometry
3099+
Point3D* p1 = mPointVector[currentEdge->getPoint()];
3100+
Point3D* p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
3101+
QgsPolyline lineGeom;
3102+
lineGeom.push_back( QgsPoint( p1->getX(), p1->getY() ) );
3103+
lineGeom.push_back( QgsPoint( p2->getX(), p2->getY() ) );
3104+
QgsGeometry* geom = QgsGeometry::fromPolyline( lineGeom );
3105+
edgeLineFeature.setGeometry( geom );
3106+
3107+
//attributes
3108+
QString attributeString;
3109+
if ( currentEdge->getForced() )
3110+
{
3111+
if ( currentEdge->getBreak() )
3112+
{
3113+
attributeString = "break line";
3114+
}
3115+
else
3116+
{
3117+
attributeString = "structure line";
3118+
}
3119+
}
3120+
edgeLineFeature.addAttribute( 0, attributeString );
3121+
3122+
writer.addFeature( edgeLineFeature );
3123+
}
3124+
alreadyVisitedEdges[i] = true;
3125+
}
3126+
3127+
return true;
3128+
}
3129+
30733130
double DualEdgeTriangulation::swapMinAngle( int edge ) const
30743131
{
30753132
Point3D* p1 = getPoint( mHalfEdge[edge]->getPoint() );

‎src/plugins/interpolation/DualEdgeTriangulation.h

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ class DualEdgeTriangulation: public Triangulation
4343
DualEdgeTriangulation();
4444
DualEdgeTriangulation( int nop, Triangulation* decorator );
4545
virtual ~DualEdgeTriangulation();
46-
/**Adds a line (e.g. a break-, structure- or an isoline) to the triangulation*/
46+
/**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*/
4747
void addLine( Line3D* line, bool breakline );
4848
/**Adds a point to the triangulation and returns the number of this point in case of success or -100 in case of failure*/
4949
int addPoint( Point3D* p );
@@ -104,6 +104,9 @@ class DualEdgeTriangulation: public Triangulation
104104
virtual bool swapEdge( double x, double y );
105105
/**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*/
106106
virtual QList<int>* getPointsAroundEdge( double x, double y );
107+
/**Saves the triangulation as a (line) shapefile
108+
@return true in case of success*/
109+
virtual bool saveAsShapefile( const QString& fileName ) const;
107110

108111
protected:
109112
/**X-coordinate of the upper right corner of the bounding box*/
@@ -137,7 +140,7 @@ class DualEdgeTriangulation: public Triangulation
137140
/**inserts an edge and makes sure, everything is ok with the storage of the edge. The number of the HalfEdge is returned*/
138141
unsigned int insertEdge( int dual, int next, int point, bool mbreak, bool forced );
139142
/**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*/
140-
//int insertForcedSegment(int p1, int p2, bool breakline);
143+
int insertForcedSegment( int p1, int p2, bool breakline );
141144
/**Treshold for the leftOfTest to handle numerical instabilities*/
142145
//const static double leftOfTresh=0.00001;
143146
/**Security to prevent endless loops in 'baseEdgeOfTriangle'. It there are more iteration then this number, the point will not be inserted*/
@@ -165,7 +168,7 @@ class DualEdgeTriangulation: public Triangulation
165168
/**Returns true, if it is possible to swap an edge, otherwise false(concave quad or edge on (or outside) the convex hull)*/
166169
bool swapPossible( unsigned int edge );
167170
/**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*/
168-
//void triangulatePolygon(QList<int>* poly, QList<int>* free, int mainedge);
171+
void triangulatePolygon( QList<int>* poly, QList<int>* free, int mainedge );
169172
/**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*/
170173
bool halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const;
171174
/**Calculates the minimum angle, which would be present, if the specified halfedge would be swapped*/

‎src/plugins/interpolation/Triangulation.h

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@ class Triangulation
3131
/**Enumeration describing the behaviour, if two forced lines cross. SnappingType_VERTICE means, that the second inserted forced line is snapped to the closest vertice of the first inserted forced line. DELETE_FIRST means, that the status of the first inserted forced line is reset to that of a normal edge (so that the second inserted forced line remain and the first not*/
3232
enum forcedCrossBehaviour {SnappingType_VERTICE, DELETE_FIRST, INSERT_VERTICE};
3333
virtual ~Triangulation();
34-
/**Adds a line (e.g. a break-, structure- or an isoline) to the triangulation*/
34+
/**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*/
3535
virtual void addLine( Line3D* line, bool breakline ) = 0;
3636
/**Adds a point to the triangulation*/
3737
virtual int addPoint( Point3D* p ) = 0;
@@ -88,6 +88,9 @@ class Triangulation
8888
//virtual bool saveToTAFF(QString fileName) const=0;
8989
/**Swaps the edge which is closest to the point with x and y coordinates (if this is possible)*/
9090
virtual bool swapEdge( double x, double y ) = 0;
91+
/**Saves the triangulation as a (line) shapefile
92+
@return true in case of success*/
93+
virtual bool saveAsShapefile( const QString& fileName ) const = 0;
9194
};
9295

9396
inline Triangulation::~Triangulation()

‎src/plugins/interpolation/qgsidwinterpolator.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,12 +19,12 @@
1919
#include <cmath>
2020
#include <limits>
2121

22-
QgsIDWInterpolator::QgsIDWInterpolator( const QList<QgsVectorLayer*>& vlayers ): QgsInterpolator( vlayers ), mDistanceCoefficient( 2.0 )
22+
QgsIDWInterpolator::QgsIDWInterpolator( const QList<LayerData>& layerData ): QgsInterpolator( layerData ), mDistanceCoefficient( 2.0 )
2323
{
2424

2525
}
2626

27-
QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList<QgsVectorLayer*>() ), mDistanceCoefficient( 2.0 )
27+
QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList<LayerData>() ), mDistanceCoefficient( 2.0 )
2828
{
2929

3030
}

‎src/plugins/interpolation/qgsidwinterpolator.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
class QgsIDWInterpolator: public QgsInterpolator
2424
{
2525
public:
26-
QgsIDWInterpolator( const QList<QgsVectorLayer*>& vlayers );
26+
QgsIDWInterpolator( const QList<LayerData>& layerData );
2727
~QgsIDWInterpolator();
2828

2929
/**Calculates interpolation value for map coordinates x, y

‎src/plugins/interpolation/qgsidwinterpolatordialog.cpp

Lines changed: 1 addition & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -13,15 +13,7 @@ QgsIDWInterpolatorDialog::~QgsIDWInterpolatorDialog()
1313

1414
QgsInterpolator* QgsIDWInterpolatorDialog::createInterpolator() const
1515
{
16-
QList<QgsVectorLayer*> inputLayerList;
17-
18-
QList< QPair <QgsVectorLayer*, QgsInterpolator::InputType> >::const_iterator data_it = mInputData.constBegin();
19-
for ( ; data_it != mInputData.constEnd(); ++data_it )
20-
{
21-
inputLayerList.push_back( data_it->first );
22-
}
23-
24-
QgsIDWInterpolator* theInterpolator = new QgsIDWInterpolator( inputLayerList );
16+
QgsIDWInterpolator* theInterpolator = new QgsIDWInterpolator( mInputData );
2517
theInterpolator->setDistanceCoefficient( mPSpinBox->value() );
2618
return theInterpolator;
2719
}

0 commit comments

Comments
 (0)
Please sign in to comment.