Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
FEATURE: possibility to add line layers as constrains for triangulati…
…on in interpolation plugin

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@11211 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mhugent committed Jul 30, 2009
1 parent 664dfad commit 83c8bd1
Show file tree
Hide file tree
Showing 16 changed files with 500 additions and 182 deletions.
127 changes: 92 additions & 35 deletions src/plugins/interpolation/DualEdgeTriangulation.cc
Expand Up @@ -17,7 +17,9 @@

#include "DualEdgeTriangulation.h"
#include <map>
#include "qgsgeometry.h"
#include "qgslogger.h"
#include "qgsvectorfilewriter.h"

double leftOfTresh = 0.00000001;

Expand Down Expand Up @@ -75,7 +77,7 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline )
for ( i = 0;i < line->getSize();i++ )
{
line->goToNext();
actpoint = mDecorator->addPoint( line->getPoint() );
actpoint = /*mDecorator->*/addPoint( line->getPoint() );
if ( actpoint != -100 )
{
i++;
Expand All @@ -85,20 +87,22 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline )

if ( actpoint == -100 )//no point of the line could be inserted
{
delete line;
return;
}

for ( ;i < line->getSize();i++ )
{
line->goToNext();
currentpoint = mDecorator->addPoint( line->getPoint() );
//if(currentpoint!=-100&&actpoint!=-100&&currentpoint!=actpoint)//-100 is the return value if the point could not be not inserted
//{
// insertForcedSegment(actpoint,currentpoint,breakline);
//}
currentpoint = /*mDecorator->*/addPoint( line->getPoint() );
if ( currentpoint != -100 && actpoint != -100 && currentpoint != actpoint )//-100 is the return value if the point could not be not inserted
{
insertForcedSegment( actpoint, currentpoint, breakline );
}
actpoint = currentpoint;
}
}
delete line;
}

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

}

#if 0
int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
{
if ( p1 == p2 )
Expand Down Expand Up @@ -1436,8 +1439,8 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
}

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

//finish the polygon on the left side
int actpointl = p2;
QListIterator<int> leftiter;
leftiter = crossedEdges.fromLast();
QList<int>::const_iterator leftiter; //todo: is there a better way to set an iterator to the last list element?
leftiter = crossedEdges.constEnd();
--leftiter;
while ( true )
{
int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext()]->getPoint();
Expand All @@ -1479,7 +1483,7 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
int theedge = mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext();
leftPolygon.append( theedge );
}
if ( leftiter == crossedEdges.begin() )
if ( leftiter == crossedEdges.constBegin() )
{break;}
--leftiter;
}
Expand All @@ -1488,9 +1492,9 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
leftPolygon.append( mHalfEdge[crossedEdges.first()]->getNext() );

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

//set the necessary nexts of leftPolygon(exept the first)
int actedgel = leftPolygon[1];
for ( leftiter = leftPolygon.at( 2 );leftiter != leftPolygon.end();++leftiter )
leftiter = leftPolygon.constBegin(); leftiter += 2;
for ( ;leftiter != leftPolygon.constEnd();++leftiter )
{
mHalfEdge[actedgel]->setNext(( *leftiter ) );
actedgel = ( *leftiter );
}

//set all the necessary nexts of rightPolygon
int actedger = rightPolygon[1];
for ( rightiter = rightPolygon.at( 2 );rightiter != rightPolygon.end();++rightiter )
rightiter = rightPolygon.constBegin(); rightiter += 2;
for ( ;rightiter != rightPolygon.constEnd();++rightiter )
{
mHalfEdge[actedger]->setNext(( *rightiter ) );
actedger = ( *( rightiter ) );
Expand All @@ -1532,8 +1538,6 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
mHalfEdge[rightPolygon.first()]->setPoint( p1 );
mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );



triangulatePolygon( &leftPolygon, &freelist, firstedge );
triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );

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

return leftPolygon.first();
}
#endif //0

void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b )
{
Expand Down Expand Up @@ -2442,7 +2445,6 @@ bool DualEdgeTriangulation::swapPossible( unsigned int edge )
return true;
}

#if 0
void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* free, int mainedge )
{
if ( poly && free )
Expand All @@ -2453,13 +2455,13 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
}

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

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

mHalfEdge[inserta]->setNext(( *( poly->at( 1 ) ) ) );
mHalfEdge[inserta]->setNext(( poly->at( 1 ) ) );
mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
mHalfEdge[insertb]->setNext( nextdistedge );
mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
mHalfEdge[distedge]->setNext( inserta );
mHalfEdge[mainedge]->setNext( insertb );

QValueList<int> polya;
for ( iterator = ( ++( poly->begin() ) );( *iterator ) != nextdistedge;++iterator )
QList<int> polya;
for ( iterator = ( ++( poly->constBegin() ) );( *iterator ) != nextdistedge;++iterator )
{
polya.append(( *iterator ) );
}
Expand All @@ -2507,16 +2509,16 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
int insertb = mHalfEdge[inserta]->getDual();
free->pop_front();

mHalfEdge[inserta]->setNext(( *( poly->at( 2 ) ) ) );
mHalfEdge[inserta]->setNext(( poly->at( 2 ) ) );
mHalfEdge[inserta]->setPoint( mHalfEdge[distedge]->getPoint() );
mHalfEdge[insertb]->setNext( mainedge );
mHalfEdge[insertb]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
mHalfEdge[distedge]->setNext( insertb );
mHalfEdge[( *( --poly->end() ) )]->setNext( inserta );

QValueList<int> polya;
iterator = poly->at( 2 );
while ( iterator != poly->end() )
QList<int> polya;
iterator = poly->constBegin(); iterator += 2;
while ( iterator != poly->constEnd() )
{
polya.append(( *iterator ) );
++iterator;
Expand All @@ -2536,7 +2538,7 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
int insertd = mHalfEdge[insertc]->getDual();
free->pop_front();

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

//build two new polygons for recursive triangulation
QValueList<int> polya;
QValueList<int> polyb;
QList<int> polya;
QList<int> polyb;

for ( iterator = ++( poly->begin() );( *iterator ) != nextdistedge;++iterator )
for ( iterator = ++( poly->constBegin() );( *iterator ) != nextdistedge;++iterator )
{
polya.append(( *iterator ) );
}
polya.prepend( inserta );


while ( iterator != poly->end() )
while ( iterator != poly->constEnd() )
{
polyb.append(( *iterator ) );
++iterator;
Expand All @@ -2578,7 +2580,6 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
}

}
#endif //0

bool DualEdgeTriangulation::pointInside( double x, double y )
{
Expand Down Expand Up @@ -3070,6 +3071,62 @@ QList<int>* DualEdgeTriangulation::getPointsAroundEdge( double x, double y )
}
}

bool DualEdgeTriangulation::saveAsShapefile( const QString& fileName ) const
{
QgsFieldMap fields;
fields.insert( 0, QgsField( "type", QVariant::String, "String" ) );
QgsVectorFileWriter writer( fileName, "Utf-8", fields, QGis::WKBLineString, 0 );
if ( writer.hasError() != QgsVectorFileWriter::NoError )
{
return false;
}

bool alreadyVisitedEdges[mHalfEdge.size()];

for ( int i = 0; i < mHalfEdge.size(); ++i )
{
alreadyVisitedEdges[i] = false;
}

for ( int i = 0; i < mHalfEdge.size(); ++i )
{
HalfEdge* currentEdge = mHalfEdge[i];
if ( currentEdge->getPoint() != -1 && mHalfEdge[currentEdge->getDual()]->getPoint() != -1 && !alreadyVisitedEdges[currentEdge->getDual()] )
{
QgsFeature edgeLineFeature;

//geometry
Point3D* p1 = mPointVector[currentEdge->getPoint()];
Point3D* p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
QgsPolyline lineGeom;
lineGeom.push_back( QgsPoint( p1->getX(), p1->getY() ) );
lineGeom.push_back( QgsPoint( p2->getX(), p2->getY() ) );
QgsGeometry* geom = QgsGeometry::fromPolyline( lineGeom );
edgeLineFeature.setGeometry( geom );

//attributes
QString attributeString;
if ( currentEdge->getForced() )
{
if ( currentEdge->getBreak() )
{
attributeString = "break line";
}
else
{
attributeString = "structure line";
}
}
edgeLineFeature.addAttribute( 0, attributeString );

writer.addFeature( edgeLineFeature );
}
alreadyVisitedEdges[i] = true;
}

return true;
}

double DualEdgeTriangulation::swapMinAngle( int edge ) const
{
Point3D* p1 = getPoint( mHalfEdge[edge]->getPoint() );
Expand Down
9 changes: 6 additions & 3 deletions src/plugins/interpolation/DualEdgeTriangulation.h
Expand Up @@ -43,7 +43,7 @@ class DualEdgeTriangulation: public Triangulation
DualEdgeTriangulation();
DualEdgeTriangulation( int nop, Triangulation* decorator );
virtual ~DualEdgeTriangulation();
/**Adds a line (e.g. a break-, structure- or an isoline) to the triangulation*/
/**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 );
/**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 );
Expand Down Expand Up @@ -104,6 +104,9 @@ class DualEdgeTriangulation: public Triangulation
virtual bool swapEdge( double x, double y );
/**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 );
/**Saves the triangulation as a (line) shapefile
@return true in case of success*/
virtual bool saveAsShapefile( const QString& fileName ) const;

protected:
/**X-coordinate of the upper right corner of the bounding box*/
Expand Down Expand Up @@ -137,7 +140,7 @@ class DualEdgeTriangulation: public Triangulation
/**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);
int insertForcedSegment( int p1, int p2, bool breakline );
/**Treshold 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*/
Expand Down Expand Up @@ -165,7 +168,7 @@ class DualEdgeTriangulation: public Triangulation
/**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);
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*/
Expand Down
5 changes: 4 additions & 1 deletion src/plugins/interpolation/Triangulation.h
Expand Up @@ -31,7 +31,7 @@ class Triangulation
/**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*/
enum forcedCrossBehaviour {SnappingType_VERTICE, DELETE_FIRST, INSERT_VERTICE};
virtual ~Triangulation();
/**Adds a line (e.g. a break-, structure- or an isoline) to the triangulation*/
/**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*/
virtual void addLine( Line3D* line, bool breakline ) = 0;
/**Adds a point to the triangulation*/
virtual int addPoint( Point3D* p ) = 0;
Expand Down Expand Up @@ -88,6 +88,9 @@ class Triangulation
//virtual bool saveToTAFF(QString fileName) const=0;
/**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 ) = 0;
/**Saves the triangulation as a (line) shapefile
@return true in case of success*/
virtual bool saveAsShapefile( const QString& fileName ) const = 0;
};

inline Triangulation::~Triangulation()
Expand Down
4 changes: 2 additions & 2 deletions src/plugins/interpolation/qgsidwinterpolator.cpp
Expand Up @@ -19,12 +19,12 @@
#include <cmath>
#include <limits>

QgsIDWInterpolator::QgsIDWInterpolator( const QList<QgsVectorLayer*>& vlayers ): QgsInterpolator( vlayers ), mDistanceCoefficient( 2.0 )
QgsIDWInterpolator::QgsIDWInterpolator( const QList<LayerData>& layerData ): QgsInterpolator( layerData ), mDistanceCoefficient( 2.0 )
{

}

QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList<QgsVectorLayer*>() ), mDistanceCoefficient( 2.0 )
QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList<LayerData>() ), mDistanceCoefficient( 2.0 )
{

}
Expand Down
2 changes: 1 addition & 1 deletion src/plugins/interpolation/qgsidwinterpolator.h
Expand Up @@ -23,7 +23,7 @@
class QgsIDWInterpolator: public QgsInterpolator
{
public:
QgsIDWInterpolator( const QList<QgsVectorLayer*>& vlayers );
QgsIDWInterpolator( const QList<LayerData>& layerData );
~QgsIDWInterpolator();

/**Calculates interpolation value for map coordinates x, y
Expand Down
10 changes: 1 addition & 9 deletions src/plugins/interpolation/qgsidwinterpolatordialog.cpp
Expand Up @@ -13,15 +13,7 @@ QgsIDWInterpolatorDialog::~QgsIDWInterpolatorDialog()

QgsInterpolator* QgsIDWInterpolatorDialog::createInterpolator() const
{
QList<QgsVectorLayer*> inputLayerList;

QList< QPair <QgsVectorLayer*, QgsInterpolator::InputType> >::const_iterator data_it = mInputData.constBegin();
for ( ; data_it != mInputData.constEnd(); ++data_it )
{
inputLayerList.push_back( data_it->first );
}

QgsIDWInterpolator* theInterpolator = new QgsIDWInterpolator( inputLayerList );
QgsIDWInterpolator* theInterpolator = new QgsIDWInterpolator( mInputData );
theInterpolator->setDistanceCoefficient( mPSpinBox->value() );
return theInterpolator;
}

0 comments on commit 83c8bd1

Please sign in to comment.