Skip to content

Commit

Permalink
Updated interpolation classes with the changes in trunk
Browse files Browse the repository at this point in the history
git-svn-id: http://svn.osgeo.org/qgis/trunk@11553 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mhugent committed Sep 5, 2009
1 parent 9fdde26 commit fe3b1d6
Show file tree
Hide file tree
Showing 15 changed files with 409 additions and 99 deletions.
1 change: 1 addition & 0 deletions src/analysis/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ QT4_WRAP_CPP(QGIS_ANALYSIS_MOC_SRCS ${QGIS_ANALYSIS_MOC_HDRS})
INCLUDE_DIRECTORIES(
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/../core/
${CMAKE_CURRENT_SOURCE_DIR}/../core/renderer
interpolation
${PROJ_INCLUDE_DIR}
${GEOS_INCLUDE_DIR}
Expand Down
152 changes: 119 additions & 33 deletions src/analysis/interpolation/DualEdgeTriangulation.cc
Original file line number Diff line number Diff line change
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 @@ -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);
//}
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 @@ -1121,6 +1125,10 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, Point3
int ptnr1 = mHalfEdge[edge1]->getPoint();
int ptnr2 = mHalfEdge[edge2]->getPoint();
int ptnr3 = mHalfEdge[edge3]->getPoint();
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
Expand Down Expand Up @@ -1174,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 @@ -1432,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 @@ -1463,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 @@ -1475,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 @@ -1484,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 @@ -1505,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 @@ -1528,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 @@ -1541,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 @@ -2438,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 @@ -2449,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 @@ -2472,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 @@ -2503,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 @@ -2532,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 @@ -2546,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 @@ -2574,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 @@ -3066,6 +3071,87 @@ QList<int>* DualEdgeTriangulation::getPointsAroundEdge( double x, double y )
}
}

bool DualEdgeTriangulation::saveAsShapefile( const QString& fileName ) const
{
QString shapeFileName = fileName;

QgsFieldMap fields;
fields.insert( 0, QgsField( "type", QVariant::String, "String" ) );

// add the extension if not present
if ( shapeFileName.indexOf( ".shp" ) == -1 )
{
shapeFileName += ".shp";
}

//delete already existing files
if ( QFile::exists( shapeFileName ) )
{
if ( !QgsVectorFileWriter::deleteShapeFile( shapeFileName ) )
{
return false;
}
}

QgsVectorFileWriter writer( shapeFileName, "Utf-8", fields, QGis::WKBLineString, 0 );
if ( writer.hasError() != QgsVectorFileWriter::NoError )
{
return false;
}

bool *alreadyVisitedEdges = new bool[mHalfEdge.size()];
if ( !alreadyVisitedEdges )
{
QgsDebugMsg( "out of memory" );
return false;
}

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;
}

delete [] alreadyVisitedEdges;

return true;
}

double DualEdgeTriangulation::swapMinAngle( int edge ) const
{
Point3D* p1 = getPoint( mHalfEdge[edge]->getPoint() );
Expand Down
14 changes: 11 additions & 3 deletions src/analysis/interpolation/DualEdgeTriangulation.h
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@ class ANALYSIS_EXPORT 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*/
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 );
/**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 +105,9 @@ class ANALYSIS_EXPORT 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 +141,7 @@ class ANALYSIS_EXPORT 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 +169,7 @@ class ANALYSIS_EXPORT 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 All @@ -188,6 +192,10 @@ inline DualEdgeTriangulation::DualEdgeTriangulation( int nop, Triangulation* dec
{
mPointVector.reserve( nop );
mHalfEdge.reserve( nop );
if ( !mDecorator )
{
mDecorator = this;
}
}

inline double DualEdgeTriangulation::getXMax() const
Expand Down
Loading

0 comments on commit fe3b1d6

Please sign in to comment.