diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index ce32dcc8b75e..e4cd86378bef 100644 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -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} diff --git a/src/analysis/interpolation/DualEdgeTriangulation.cc b/src/analysis/interpolation/DualEdgeTriangulation.cc index 3bfc8b91b82f..2d8ba2cc902a 100644 --- a/src/analysis/interpolation/DualEdgeTriangulation.cc +++ b/src/analysis/interpolation/DualEdgeTriangulation.cc @@ -17,7 +17,9 @@ #include "DualEdgeTriangulation.h" #include +#include "qgsgeometry.h" #include "qgslogger.h" +#include "qgsvectorfilewriter.h" double leftOfTresh = 0.00000001; @@ -85,6 +87,7 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline ) if ( actpoint == -100 )//no point of the line could be inserted { + delete line; return; } @@ -92,13 +95,14 @@ void DualEdgeTriangulation::addLine( Line3D* line, bool breakline ) { line->goToNext(); currentpoint = mDecorator->addPoint( line->getPoint() ); - //if(currentpoint!=-100&&actpoint!=-100&¤tpoint!=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 ) @@ -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() ); @@ -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 ) @@ -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 iter; - for ( iter = crossedEdges.begin();iter != crossedEdges.end();++iter ) + QList::const_iterator iter; + for ( iter = crossedEdges.constBegin();iter != crossedEdges.constEnd();++iter ) { mHalfEdge[( *( iter ) )]->setForced( false ); mHalfEdge[( *( iter ) )]->setBreak( false ); @@ -1463,8 +1470,9 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline ) //finish the polygon on the left side int actpointl = p2; - QListIterator leftiter; - leftiter = crossedEdges.fromLast(); + QList::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(); @@ -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; } @@ -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 rightiter; + QList::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 ) @@ -1505,7 +1513,8 @@ 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 ); @@ -1513,7 +1522,8 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline ) //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 ) ); @@ -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 ); @@ -1541,7 +1549,6 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline ) return leftPolygon.first(); } -#endif //0 void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b ) { @@ -2438,7 +2445,6 @@ bool DualEdgeTriangulation::swapPossible( unsigned int edge ) return true; } -#if 0 void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* free, int mainedge ) { if ( poly && free ) @@ -2449,13 +2455,13 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* fr } //search for the edge pointing on the closest point(distedge) and for the next(nextdistedge) - QValueListIterator iterator = ++( poly->begin() );//go to the second edge + QList::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 ) { @@ -2472,15 +2478,15 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* 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 polya; - for ( iterator = ( ++( poly->begin() ) );( *iterator ) != nextdistedge;++iterator ) + QList polya; + for ( iterator = ( ++( poly->constBegin() ) );( *iterator ) != nextdistedge;++iterator ) { polya.append(( *iterator ) ); } @@ -2503,16 +2509,16 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* 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 polya; - iterator = poly->at( 2 ); - while ( iterator != poly->end() ) + QList polya; + iterator = poly->constBegin(); iterator += 2; + while ( iterator != poly->constEnd() ) { polya.append(( *iterator ) ); ++iterator; @@ -2532,7 +2538,7 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* 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() ); @@ -2546,17 +2552,17 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* fr mHalfEdge[( *( --poly->end() ) )]->setNext( insertc ); //build two new polygons for recursive triangulation - QValueList polya; - QValueList polyb; + QList polya; + QList 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; @@ -2574,7 +2580,6 @@ void DualEdgeTriangulation::triangulatePolygon( QList* poly, QList* fr } } -#endif //0 bool DualEdgeTriangulation::pointInside( double x, double y ) { @@ -3066,6 +3071,87 @@ QList* 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() ); diff --git a/src/analysis/interpolation/DualEdgeTriangulation.h b/src/analysis/interpolation/DualEdgeTriangulation.h index b15036f8d121..7d2e895448fd 100644 --- a/src/analysis/interpolation/DualEdgeTriangulation.h +++ b/src/analysis/interpolation/DualEdgeTriangulation.h @@ -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 ); @@ -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* 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 +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*/ @@ -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* poly, QList* free, int mainedge); + void triangulatePolygon( QList* poly, QList* 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*/ @@ -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 diff --git a/src/analysis/interpolation/MathUtils.h b/src/analysis/interpolation/MathUtils.h index 6a35207de423..56da8dd104d2 100644 --- a/src/analysis/interpolation/MathUtils.h +++ b/src/analysis/interpolation/MathUtils.h @@ -17,12 +17,16 @@ #ifndef MATHUTILS_H #define MATHUTILS_H +#ifndef Q_OS_MACX #include +#else +#include +#endif #include "Vector3D.h" #include "Point3D.h" -namespace MathUtils +namespace ANALYSIS_EXPORT MathUtils { /**Calculates the barycentric coordinates of a point (x,y) with respect to p1, p2, p3 and stores the three barycentric coordinates in 'result'. Thus the u-coordinate is stored in result::x, the v-coordinate in result::y and the w-coordinate in result::z. Attention: p1, p2 and p3 have to be ordered counterclockwise*/ bool calcBarycentricCoordinates( double x, double y, Point3D* p1, Point3D* p2, Point3D* p3, Point3D* result ); diff --git a/src/analysis/interpolation/Point3D.h b/src/analysis/interpolation/Point3D.h index 34abc895f993..3c069fc4e595 100644 --- a/src/analysis/interpolation/Point3D.h +++ b/src/analysis/interpolation/Point3D.h @@ -17,7 +17,11 @@ #ifndef POINT3D_H #define POINT3D_H +#ifndef Q_OS_MACX #include +#else +#include +#endif #include /**Point3D is a class to represent a three dimensional point*/ diff --git a/src/analysis/interpolation/Triangulation.h b/src/analysis/interpolation/Triangulation.h index 5d25c04d706c..680f027335cd 100644 --- a/src/analysis/interpolation/Triangulation.h +++ b/src/analysis/interpolation/Triangulation.h @@ -31,7 +31,7 @@ class ANALYSIS_EXPORT 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 ANALYSIS_EXPORT 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() diff --git a/src/analysis/interpolation/Vector3D.h b/src/analysis/interpolation/Vector3D.h index 6a19652042bd..d67df75cbd77 100644 --- a/src/analysis/interpolation/Vector3D.h +++ b/src/analysis/interpolation/Vector3D.h @@ -17,12 +17,16 @@ #ifndef VECTOR3D_H #define VECTOR3D_H +#ifndef Q_OS_MACX #include +#else +#include +#endif class ANALYSIS_EXPORT Vector3D /** - Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values. In fact, the class is the same as Point3D. The name 'vector' makes it easier to understand the programms. - */ + Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values. In fact, the class is the same as Point3D. The name 'vector' makes it easier to understand the programms. + */ { protected: diff --git a/src/analysis/interpolation/qgsgridfilewriter.cpp b/src/analysis/interpolation/qgsgridfilewriter.cpp index ee67eed7a392..a5e9be713656 100644 --- a/src/analysis/interpolation/qgsgridfilewriter.cpp +++ b/src/analysis/interpolation/qgsgridfilewriter.cpp @@ -20,10 +20,10 @@ #include #include -QgsGridFileWriter::QgsGridFileWriter( QgsInterpolator* i, QString outputPath, QgsRectangle extent, int nCols, int nRows ): mInterpolator( i ), mOutputFilePath( outputPath ), mInterpolationExtent( extent ), mNumColumns( nCols ), mNumRows( nRows ) +QgsGridFileWriter::QgsGridFileWriter( QgsInterpolator* i, QString outputPath, QgsRectangle extent, int nCols, int nRows , double cellSizeX, double cellSizeY ): \ + mInterpolator( i ), mOutputFilePath( outputPath ), mInterpolationExtent( extent ), mNumColumns( nCols ), mNumRows( nRows ), mCellSizeX( cellSizeX ), mCellSizeY( cellSizeY ) { - mCellSizeX = ( mInterpolationExtent.xMaximum() - mInterpolationExtent.xMinimum() ) / mNumColumns; - mCellSizeY = ( mInterpolationExtent.yMaximum() - mInterpolationExtent.yMinimum() ) / mNumRows; + } QgsGridFileWriter::QgsGridFileWriter(): mInterpolator( 0 ) @@ -65,6 +65,7 @@ int QgsGridFileWriter::writeFile( bool showProgressDialog ) progressDialog = new QProgressDialog( QObject::tr( "Interpolating..." ), QObject::tr( "Abort" ), 0, mNumRows, 0 ); progressDialog->setWindowModality( Qt::WindowModal ); } + for ( int i = 0; i < mNumRows; ++i ) { currentXValue = mInterpolationExtent.xMinimum(); diff --git a/src/analysis/interpolation/qgsgridfilewriter.h b/src/analysis/interpolation/qgsgridfilewriter.h index 6c440f49404e..602f7b73eae0 100644 --- a/src/analysis/interpolation/qgsgridfilewriter.h +++ b/src/analysis/interpolation/qgsgridfilewriter.h @@ -29,7 +29,7 @@ class QgsInterpolator; class ANALYSIS_EXPORT QgsGridFileWriter { public: - QgsGridFileWriter( QgsInterpolator* i, QString outputPath, QgsRectangle extent, int nCols, int nRows ); + QgsGridFileWriter( QgsInterpolator* i, QString outputPath, QgsRectangle extent, int nCols, int nRows, double cellSizeX, double cellSizeY ); ~QgsGridFileWriter(); /**Writes the grid file. diff --git a/src/analysis/interpolation/qgsidwinterpolator.cpp b/src/analysis/interpolation/qgsidwinterpolator.cpp index 19910d9ab4f5..d75b4584da52 100644 --- a/src/analysis/interpolation/qgsidwinterpolator.cpp +++ b/src/analysis/interpolation/qgsidwinterpolator.cpp @@ -19,12 +19,12 @@ #include #include -QgsIDWInterpolator::QgsIDWInterpolator( const QList& vlayers ): QgsInterpolator( vlayers ), mDistanceCoefficient( 2.0 ) +QgsIDWInterpolator::QgsIDWInterpolator( const QList& layerData ): QgsInterpolator( layerData ), mDistanceCoefficient( 2.0 ) { } -QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList() ), mDistanceCoefficient( 2.0 ) +QgsIDWInterpolator::QgsIDWInterpolator(): QgsInterpolator( QList() ), mDistanceCoefficient( 2.0 ) { } @@ -62,6 +62,11 @@ int QgsIDWInterpolator::interpolatePoint( double x, double y, double& result ) sumDenominator += currentWeight; } + if ( sumDenominator == 0.0 ) + { + return 1; + } + result = sumCounter / sumDenominator; return 0; } diff --git a/src/analysis/interpolation/qgsidwinterpolator.h b/src/analysis/interpolation/qgsidwinterpolator.h index 3fe87c5be2b1..68d2d0713e97 100644 --- a/src/analysis/interpolation/qgsidwinterpolator.h +++ b/src/analysis/interpolation/qgsidwinterpolator.h @@ -23,7 +23,7 @@ class ANALYSIS_EXPORT QgsIDWInterpolator: public QgsInterpolator { public: - QgsIDWInterpolator( const QList& vlayers ); + QgsIDWInterpolator( const QList& layerData ); ~QgsIDWInterpolator(); /**Calculates interpolation value for map coordinates x, y diff --git a/src/analysis/interpolation/qgsinterpolator.cpp b/src/analysis/interpolation/qgsinterpolator.cpp index ac53d35a66f4..56e2be6dc3d2 100644 --- a/src/analysis/interpolation/qgsinterpolator.cpp +++ b/src/analysis/interpolation/qgsinterpolator.cpp @@ -23,12 +23,12 @@ #else #include #endif -#ifdef _MSC_VER +#ifdef WIN32 #include #define isnan(f) _isnan(f) #endif -QgsInterpolator::QgsInterpolator( const QList& vlayers ): mDataIsCached( false ), mVectorLayers( vlayers ), zCoordInterpolation( false ), mValueAttribute( -1 ) +QgsInterpolator::QgsInterpolator( const QList& layerData ): mDataIsCached( false ), mLayerData( layerData ) { } @@ -43,15 +43,9 @@ QgsInterpolator::~QgsInterpolator() } -void QgsInterpolator::enableAttributeValueInterpolation( int attribute ) -{ - mValueAttribute = attribute; - zCoordInterpolation = false; -} - int QgsInterpolator::cacheBaseData() { - if ( mVectorLayers.size() < 1 ) + if ( mLayerData.size() < 1 ) { return 0; } @@ -60,25 +54,25 @@ int QgsInterpolator::cacheBaseData() mCachedBaseData.clear(); mCachedBaseData.reserve( 100000 ); - QList::iterator v_it = mVectorLayers.begin(); + QList::iterator v_it = mLayerData.begin(); - for ( ; v_it != mVectorLayers.end(); ++v_it ) + for ( ; v_it != mLayerData.end(); ++v_it ) { - if (( *v_it ) == 0 ) + if ( v_it->vectorLayer == 0 ) { continue; } - QgsVectorDataProvider* provider = ( *v_it )->dataProvider(); + QgsVectorDataProvider* provider = v_it->vectorLayer->dataProvider(); if ( !provider ) { return 2; } QgsAttributeList attList; - if ( !zCoordInterpolation ) + if ( !v_it->zCoordInterpolation ) { - attList.push_back( mValueAttribute ); + attList.push_back( v_it->interpolationAttribute ); } provider->select( attList ); @@ -89,22 +83,22 @@ int QgsInterpolator::cacheBaseData() while ( provider->nextFeature( theFeature ) ) { - if ( !zCoordInterpolation ) + if ( !v_it->zCoordInterpolation ) { QgsAttributeMap attMap = theFeature.attributeMap(); - QgsAttributeMap::const_iterator att_it = attMap.find( mValueAttribute ); - if ( att_it == attMap.end() ) //attribute not found, something must be wrong + QgsAttributeMap::const_iterator att_it = attMap.find( v_it->interpolationAttribute ); + if ( att_it == attMap.end() ) //attribute not found, something must be wrong (e.g. NULL value) { - return 3; + continue; } - attributeValue = att_it.value().toDouble(&attributeConversionOk); - if(!attributeConversionOk || isnan(attributeValue)) //don't consider vertices with attributes like 'nan' for the interpolation + attributeValue = att_it.value().toDouble( &attributeConversionOk ); + if ( !attributeConversionOk || isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation { continue; } } - if ( addVerticesToCache( theFeature.geometry(), attributeValue ) != 0 ) + if ( addVerticesToCache( theFeature.geometry(), v_it->zCoordInterpolation, attributeValue ) != 0 ) { return 3; } @@ -114,7 +108,7 @@ int QgsInterpolator::cacheBaseData() return 0; } -int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, double attributeValue ) +int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, bool zCoord, double attributeValue ) { if ( !geom ) { @@ -136,7 +130,7 @@ int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, double attributeValu theVertex.x = *(( double * )( currentWkbPtr ) ); currentWkbPtr += sizeof( double ); theVertex.y = *(( double * )( currentWkbPtr ) ); - if ( zCoordInterpolation && hasZValue ) + if ( zCoord && hasZValue ) { currentWkbPtr += sizeof( double ); theVertex.z = *(( double * )( currentWkbPtr ) ); @@ -161,7 +155,7 @@ int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, double attributeValu currentWkbPtr += sizeof( double ); theVertex.y = *(( double * )( currentWkbPtr ) ); currentWkbPtr += sizeof( double ); - if ( zCoordInterpolation && hasZValue ) //skip z-coordinate for 25D geometries + if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries { theVertex.z = *(( double * )( currentWkbPtr ) ); currentWkbPtr += sizeof( double ); @@ -365,5 +359,6 @@ int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, double attributeValu default: break; } + mDataIsCached = true; return 0; } diff --git a/src/analysis/interpolation/qgsinterpolator.h b/src/analysis/interpolation/qgsinterpolator.h index a1257f4c4931..9712097c98d6 100644 --- a/src/analysis/interpolation/qgsinterpolator.h +++ b/src/analysis/interpolation/qgsinterpolator.h @@ -44,7 +44,16 @@ class ANALYSIS_EXPORT QgsInterpolator BREAK_LINES }; - QgsInterpolator( const QList& vlayers ); + /**A layer together with the information about interpolation attribute / z-coordinate interpolation and the type (point, structure line, breakline)*/ + struct LayerData + { + QgsVectorLayer* vectorLayer; + bool zCoordInterpolation; + int interpolationAttribute; + InputType mInputType; + }; + + QgsInterpolator( const QList& layerData ); virtual ~QgsInterpolator(); @@ -55,17 +64,9 @@ class ANALYSIS_EXPORT QgsInterpolator @return 0 in case of success*/ virtual int interpolatePoint( double x, double y, double& result ) = 0; - /**Use the z-coord of 25D for interpolation*/ - void enableZCoordInterpolation() {zCoordInterpolation = true;} - /**Use a vector attribute as interpolation value*/ void enableAttributeValueInterpolation( int attribute ); - /**Creates a vector layer that can be added to the main map, e.g. TIN edges for triangle interpolation. - Mouse clicks in the main map can be tracked and used for configuration (e.g. edge swaping). Default implementation does - nothing.*/ - virtual QgsVectorLayer* createVectorLayer() const {return 0;} - protected: /**Caches the vertex and value data from the provider. All the vertex data will be held in virtual memory @@ -77,18 +78,19 @@ class ANALYSIS_EXPORT QgsInterpolator /**Flag that tells if the cache already has been filled*/ bool mDataIsCached; + //Information about the input vector layers and the attributes (or z-values) that are used for interpolation + QList mLayerData; + private: QgsInterpolator(); //forbidden /**Helper method that adds the vertices of a geometry to the mCachedBaseData @param geom the geometry + @param zCoord true if the z-coordinate of the geometry is to be interpolated @param attributeValue the attribute value for interpolation (if not interpolated from z-coordinate) @return 0 in case of success*/ - int addVerticesToCache( QgsGeometry* geom, double attributeValue ); + int addVerticesToCache( QgsGeometry* geom, bool zCoord, double attributeValue ); - QList mVectorLayers; - bool zCoordInterpolation; - int mValueAttribute; }; #endif diff --git a/src/analysis/interpolation/qgstininterpolator.cpp b/src/analysis/interpolation/qgstininterpolator.cpp index 125b9c318723..3a5c682f7501 100644 --- a/src/analysis/interpolation/qgstininterpolator.cpp +++ b/src/analysis/interpolation/qgstininterpolator.cpp @@ -19,10 +19,24 @@ #include "DualEdgeTriangulation.h" #include "LinTriangleInterpolator.h" #include "Point3D.h" -//#include "qgssinglesymbolrenderer.h" +#include "qgsfeature.h" +#include "qgsgeometry.h" +#include "qgssinglesymbolrenderer.h" #include "qgsvectorlayer.h" +#include -QgsTINInterpolator::QgsTINInterpolator( const QList& inputData ): QgsInterpolator( inputData ), mTriangulation( 0 ), mTriangleInterpolator( 0 ), mIsInitialized( false ) +#ifndef Q_OS_MACX +#include +#else +#include +#endif +#ifdef WIN32 +#include +#define isnan(f) _isnan(f) +#endif + +QgsTINInterpolator::QgsTINInterpolator( const QList& inputData, bool showProgressDialog ): QgsInterpolator( inputData ), mTriangulation( 0 ), \ + mTriangleInterpolator( 0 ), mIsInitialized( false ), mShowProgressDialog( showProgressDialog ), mExportTriangulationToFile( false ) { } @@ -55,25 +69,190 @@ int QgsTINInterpolator::interpolatePoint( double x, double y, double& result ) void QgsTINInterpolator::initialize() { - if ( !mDataIsCached ) + DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, 0 ); + mTriangulation = theDualEdgeTriangulation; + + //get number of features if we use a progress bar + int nFeatures = 0; + int nProcessedFeatures = 0; + if ( mShowProgressDialog ) { - cacheBaseData(); + QList::iterator layerDataIt = mLayerData.begin(); + for ( ; layerDataIt != mLayerData.end(); ++layerDataIt ) + { + if ( layerDataIt->vectorLayer ) + { + nFeatures += layerDataIt->vectorLayer->featureCount(); + } + } } - //create DualEdgeTriangulation + QProgressDialog* theProgressDialog = 0; + if ( mShowProgressDialog ) + { + theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 ); + theProgressDialog->setWindowModality( Qt::WindowModal ); + } - DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( mCachedBaseData.size(), 0 ); - mTriangulation = theDualEdgeTriangulation; - //add all the vertices to the triangulation - QVector::const_iterator vertex_it = mCachedBaseData.constBegin(); - for ( ; vertex_it != mCachedBaseData.constEnd(); ++vertex_it ) + QgsFeature f; + QList::iterator layerDataIt = mLayerData.begin(); + for ( ; layerDataIt != mLayerData.end(); ++layerDataIt ) { - Point3D* thePoint = new Point3D( vertex_it->x, vertex_it->y, vertex_it->z ); - mTriangulation->addPoint( thePoint ); + if ( layerDataIt->vectorLayer ) + { + QgsAttributeList attList; + if ( !layerDataIt->zCoordInterpolation ) + { + attList.push_back( layerDataIt->interpolationAttribute ); + } + layerDataIt->vectorLayer->select( attList ); + while ( layerDataIt->vectorLayer->nextFeature( f ) ) + { + if ( mShowProgressDialog ) + { + if ( theProgressDialog->wasCanceled() ) + { + break; + } + theProgressDialog->setValue( nProcessedFeatures ); + } + insertData( &f, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->mInputType ); + ++nProcessedFeatures; + } + } } + delete theProgressDialog; mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation ); - mIsInitialized = true; + + //debug + if ( mExportTriangulationToFile ) + { + theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath ); + } } + +int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type ) +{ + if ( !f ) + { + return 1; + } + + QgsGeometry* g = f->geometry(); + { + if ( !g ) + { + return 2; + } + } + + //check attribute value + double attributeValue = 0; + bool attributeConversionOk = false; + if ( !zCoord ) + { + QgsAttributeMap attMap = f->attributeMap(); + QgsAttributeMap::const_iterator att_it = attMap.find( attr ); + if ( att_it == attMap.end() ) //attribute not found, something must be wrong (e.g. NULL value) + { + return 3; + } + attributeValue = att_it.value().toDouble( &attributeConversionOk ); + if ( !attributeConversionOk || isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation + { + return 4; + } + } + + //parse WKB. We cannot use the methods with QgsPoint because they don't contain z-values for 25D types + bool hasZValue = false; + double x, y, z; + unsigned char* currentWkbPtr = g->asWkb(); + + QGis::WkbType wkbType = g->wkbType(); + switch ( wkbType ) + { + case QGis::WKBPoint25D: + hasZValue = true; + case QGis::WKBPoint: + { + currentWkbPtr += ( 1 + sizeof( int ) ); + x = *(( double * )( currentWkbPtr ) ); + currentWkbPtr += sizeof( double ); + y = *(( double * )( currentWkbPtr ) ); + if ( zCoord && hasZValue ) + { + currentWkbPtr += sizeof( double ); + z = *(( double * )( currentWkbPtr ) ); + } + else + { + z = attributeValue; + } + Point3D* thePoint = new Point3D( x, y, z ); + if ( mTriangulation->addPoint( thePoint ) == -100 ) + { + return -1; + } + break; + } + case QGis::WKBLineString25D: + hasZValue = true; + case QGis::WKBLineString: + { + //maybe a structure or break line + Line3D* line = 0; + if ( type != POINTS ) + { + line = new Line3D(); + } + currentWkbPtr += ( 1 + sizeof( int ) ); + int* npoints = ( int* )currentWkbPtr; + currentWkbPtr += sizeof( int ); + for ( int index = 0; index < *npoints; ++index ) + { + x = *(( double * )( currentWkbPtr ) ); + currentWkbPtr += sizeof( double ); + y = *(( double * )( currentWkbPtr ) ); + currentWkbPtr += sizeof( double ); + if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries + { + z = *(( double * )( currentWkbPtr ) ); + } + else + { + z = attributeValue; + } + if ( hasZValue ) + { + currentWkbPtr += sizeof( double ); + } + + if ( type == POINTS ) + { + //todo: handle error code -100 + mTriangulation->addPoint( new Point3D( x, y, z ) ); + } + else + { + line->insertPoint( new Point3D( x, y, z ) ); + } + } + + if ( type != POINTS ) + { + mTriangulation->addLine( line, type == BREAK_LINES ); + } + break; + } + default: + //todo: add the same for multiline, polygon, multipolygon + break; + } + + return 0; +} + diff --git a/src/analysis/interpolation/qgstininterpolator.h b/src/analysis/interpolation/qgstininterpolator.h index 32e58f47c68c..b6817a90d8b3 100644 --- a/src/analysis/interpolation/qgstininterpolator.h +++ b/src/analysis/interpolation/qgstininterpolator.h @@ -19,15 +19,17 @@ #define QGSTININTERPOLATOR_H #include "qgsinterpolator.h" +#include class Triangulation; class TriangleInterpolator; +class QgsFeature; /**Interpolation in a triangular irregular network*/ class ANALYSIS_EXPORT QgsTINInterpolator: public QgsInterpolator { public: - QgsTINInterpolator( const QList& inputData ); + QgsTINInterpolator( const QList& inputData, bool showProgressDialog = false ); ~QgsTINInterpolator(); /**Calculates interpolation value for map coordinates x, y @@ -37,12 +39,28 @@ class ANALYSIS_EXPORT QgsTINInterpolator: public QgsInterpolator @return 0 in case of success*/ int interpolatePoint( double x, double y, double& result ); + void setExportTriangulationToFile( bool e ) {mExportTriangulationToFile = e;} + void setTriangulationFilePath( const QString& filepath ) {mTriangulationFilePath = filepath;} + private: Triangulation* mTriangulation; TriangleInterpolator* mTriangleInterpolator; bool mIsInitialized; + bool mShowProgressDialog; + /**If true: export triangulation to shapefile after initialisation*/ + bool mExportTriangulationToFile; + /**File path to export the triangulation*/ + QString mTriangulationFilePath; + /**Create dual edge triangulation*/ void initialize(); + /**Inserts the vertices of a geometry into the triangulation + @param g the geometry + @param zCoord true if the z coordinate is the interpolation attribute + @param attr interpolation attribute index (if zCoord is false) + @param type point/structure line, break line + @return 0 in case of success, -1 if the feature could not be inserted because of numerical problems*/ + int insertData( QgsFeature* f, bool zCoord, int attr, InputType type ); }; #endif