Skip to content
Permalink
Browse files

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
mhugent committed Jul 30, 2009
1 parent 664dfad commit 83c8bd1fe3d8c1cf7e53093bd037ec3d3695f685
@@ -17,7 +17,9 @@

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

double leftOfTresh = 0.00000001;

@@ -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++;
@@ -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 )
@@ -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 )
@@ -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 );
@@ -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();
@@ -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;
}
@@ -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 )
@@ -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 ) );
@@ -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 );

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

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

void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b )
{
@@ -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 )
@@ -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 )
{
@@ -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 ) );
}
@@ -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;
@@ -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() );
@@ -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;
@@ -2578,7 +2580,6 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* fr
}

}
#endif //0

bool DualEdgeTriangulation::pointInside( double x, double y )
{
@@ -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() );
@@ -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 );
@@ -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*/
@@ -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*/
@@ -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*/
@@ -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;
@@ -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()
@@ -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 )
{

}
@@ -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
@@ -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.
You can’t perform that action at this time.