[QGIS Commit] r11553 - in trunk/qgis/src/analysis: . interpolation
svn_qgis at osgeo.org
svn_qgis at osgeo.org
Sat Sep 5 05:12:26 EDT 2009
Author: mhugent
Date: 2009-09-05 05:12:26 -0400 (Sat, 05 Sep 2009)
New Revision: 11553
Modified:
trunk/qgis/src/analysis/CMakeLists.txt
trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc
trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.h
trunk/qgis/src/analysis/interpolation/MathUtils.h
trunk/qgis/src/analysis/interpolation/Point3D.h
trunk/qgis/src/analysis/interpolation/Triangulation.h
trunk/qgis/src/analysis/interpolation/Vector3D.h
trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.cpp
trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.h
trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.cpp
trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.h
trunk/qgis/src/analysis/interpolation/qgsinterpolator.cpp
trunk/qgis/src/analysis/interpolation/qgsinterpolator.h
trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp
trunk/qgis/src/analysis/interpolation/qgstininterpolator.h
Log:
Updated interpolation classes with the changes in trunk
Modified: trunk/qgis/src/analysis/CMakeLists.txt
===================================================================
--- trunk/qgis/src/analysis/CMakeLists.txt 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/CMakeLists.txt 2009-09-05 09:12:26 UTC (rev 11553)
@@ -30,6 +30,7 @@
INCLUDE_DIRECTORIES(
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/../core/
+ ${CMAKE_CURRENT_SOURCE_DIR}/../core/renderer
interpolation
${PROJ_INCLUDE_DIR}
${GEOS_INCLUDE_DIR}
Modified: trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc
===================================================================
--- trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.cc 2009-09-05 09:12:26 UTC (rev 11553)
@@ -17,7 +17,9 @@
#include "DualEdgeTriangulation.h"
#include <map>
+#include "qgsgeometry.h"
#include "qgslogger.h"
+#include "qgsvectorfilewriter.h"
double leftOfTresh = 0.00000001;
@@ -85,6 +87,7 @@
if ( actpoint == -100 )//no point of the line could be inserted
{
+ delete line;
return;
}
@@ -92,13 +95,14 @@
{
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 @@
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 @@
}
-#if 0
int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
{
if ( p1 == p2 )
@@ -1432,8 +1439,8 @@
}
//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 );
@@ -1463,8 +1470,9 @@
//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();
@@ -1475,7 +1483,7 @@
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 @@
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 )
@@ -1505,7 +1513,8 @@
//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 @@
//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 @@
mHalfEdge[rightPolygon.first()]->setPoint( p1 );
mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
-
-
triangulatePolygon( &leftPolygon, &freelist, firstedge );
triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );
@@ -1541,7 +1549,6 @@
return leftPolygon.first();
}
-#endif //0
void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b )
{
@@ -2438,7 +2445,6 @@
return true;
}
-#if 0
void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* free, int mainedge )
{
if ( poly && free )
@@ -2449,13 +2455,13 @@
}
//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 )
{
@@ -2472,15 +2478,15 @@
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 ) );
}
@@ -2503,16 +2509,16 @@
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;
@@ -2532,7 +2538,7 @@
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 @@
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;
@@ -2574,7 +2580,6 @@
}
}
-#endif //0
bool DualEdgeTriangulation::pointInside( double x, double y )
{
@@ -3066,6 +3071,87 @@
}
}
+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() );
Modified: trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/DualEdgeTriangulation.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -43,7 +43,8 @@
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 @@
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 +141,7 @@
/**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 @@
/**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*/
@@ -188,6 +192,10 @@
{
mPointVector.reserve( nop );
mHalfEdge.reserve( nop );
+ if ( !mDecorator )
+ {
+ mDecorator = this;
+ }
}
inline double DualEdgeTriangulation::getXMax() const
Modified: trunk/qgis/src/analysis/interpolation/MathUtils.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/MathUtils.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/MathUtils.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -17,12 +17,16 @@
#ifndef MATHUTILS_H
#define MATHUTILS_H
+#ifndef Q_OS_MACX
#include <cmath>
+#else
+#include <math.h>
+#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 );
Modified: trunk/qgis/src/analysis/interpolation/Point3D.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/Point3D.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/Point3D.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -17,7 +17,11 @@
#ifndef POINT3D_H
#define POINT3D_H
+#ifndef Q_OS_MACX
#include <cmath>
+#else
+#include <math.h>
+#endif
#include <iostream>
/**Point3D is a class to represent a three dimensional point*/
Modified: trunk/qgis/src/analysis/interpolation/Triangulation.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/Triangulation.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/Triangulation.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -31,7 +31,7 @@
/**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 @@
//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()
Modified: trunk/qgis/src/analysis/interpolation/Vector3D.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/Vector3D.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/Vector3D.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -17,12 +17,16 @@
#ifndef VECTOR3D_H
#define VECTOR3D_H
+#ifndef Q_OS_MACX
#include <cmath>
+#else
+#include <math.h>
+#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:
Modified: trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.cpp
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.cpp 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.cpp 2009-09-05 09:12:26 UTC (rev 11553)
@@ -20,10 +20,10 @@
#include <QFile>
#include <QProgressDialog>
-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 @@
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();
Modified: trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsgridfilewriter.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -29,7 +29,7 @@
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.
Modified: trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.cpp
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.cpp 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.cpp 2009-09-05 09:12:26 UTC (rev 11553)
@@ -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 )
{
}
@@ -62,6 +62,11 @@
sumDenominator += currentWeight;
}
+ if ( sumDenominator == 0.0 )
+ {
+ return 1;
+ }
+
result = sumCounter / sumDenominator;
return 0;
}
Modified: trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsidwinterpolator.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -23,7 +23,7 @@
class ANALYSIS_EXPORT QgsIDWInterpolator: public QgsInterpolator
{
public:
- QgsIDWInterpolator( const QList<QgsVectorLayer*>& vlayers );
+ QgsIDWInterpolator( const QList<LayerData>& layerData );
~QgsIDWInterpolator();
/**Calculates interpolation value for map coordinates x, y
Modified: trunk/qgis/src/analysis/interpolation/qgsinterpolator.cpp
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsinterpolator.cpp 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsinterpolator.cpp 2009-09-05 09:12:26 UTC (rev 11553)
@@ -23,12 +23,12 @@
#else
#include <math.h>
#endif
-#ifdef _MSC_VER
+#ifdef WIN32
#include <float.h>
#define isnan(f) _isnan(f)
#endif
-QgsInterpolator::QgsInterpolator( const QList<QgsVectorLayer*>& vlayers ): mDataIsCached( false ), mVectorLayers( vlayers ), zCoordInterpolation( false ), mValueAttribute( -1 )
+QgsInterpolator::QgsInterpolator( const QList<LayerData>& layerData ): mDataIsCached( false ), mLayerData( layerData )
{
}
@@ -43,15 +43,9 @@
}
-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 @@
mCachedBaseData.clear();
mCachedBaseData.reserve( 100000 );
- QList<QgsVectorLayer*>::iterator v_it = mVectorLayers.begin();
+ QList<LayerData>::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 @@
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 @@
return 0;
}
-int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, double attributeValue )
+int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, bool zCoord, double attributeValue )
{
if ( !geom )
{
@@ -136,7 +130,7 @@
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 @@
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 @@
default:
break;
}
+ mDataIsCached = true;
return 0;
}
Modified: trunk/qgis/src/analysis/interpolation/qgsinterpolator.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgsinterpolator.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgsinterpolator.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -44,8 +44,17 @@
BREAK_LINES
};
- QgsInterpolator( const QList<QgsVectorLayer*>& 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>& layerData );
+
virtual ~QgsInterpolator();
/**Calculates interpolation value for map coordinates x, y
@@ -55,17 +64,9 @@
@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 @@
/**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<LayerData> 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<QgsVectorLayer*> mVectorLayers;
- bool zCoordInterpolation;
- int mValueAttribute;
};
#endif
Modified: trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgstininterpolator.cpp 2009-09-05 09:12:26 UTC (rev 11553)
@@ -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 <QProgressDialog>
-QgsTINInterpolator::QgsTINInterpolator( const QList<QgsVectorLayer*>& inputData ): QgsInterpolator( inputData ), mTriangulation( 0 ), mTriangleInterpolator( 0 ), mIsInitialized( false )
+#ifndef Q_OS_MACX
+#include <cmath>
+#else
+#include <math.h>
+#endif
+#ifdef WIN32
+#include <float.h>
+#define isnan(f) _isnan(f)
+#endif
+
+QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, bool showProgressDialog ): QgsInterpolator( inputData ), mTriangulation( 0 ), \
+ mTriangleInterpolator( 0 ), mIsInitialized( false ), mShowProgressDialog( showProgressDialog ), mExportTriangulationToFile( false )
{
}
@@ -55,25 +69,190 @@
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<LayerData>::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<vertexData>::const_iterator vertex_it = mCachedBaseData.constBegin();
- for ( ; vertex_it != mCachedBaseData.constEnd(); ++vertex_it )
+ QgsFeature f;
+ QList<LayerData>::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;
- 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;
+}
+
Modified: trunk/qgis/src/analysis/interpolation/qgstininterpolator.h
===================================================================
--- trunk/qgis/src/analysis/interpolation/qgstininterpolator.h 2009-09-05 08:45:13 UTC (rev 11552)
+++ trunk/qgis/src/analysis/interpolation/qgstininterpolator.h 2009-09-05 09:12:26 UTC (rev 11553)
@@ -19,15 +19,17 @@
#define QGSTININTERPOLATOR_H
#include "qgsinterpolator.h"
+#include <QString>
class Triangulation;
class TriangleInterpolator;
+class QgsFeature;
/**Interpolation in a triangular irregular network*/
class ANALYSIS_EXPORT QgsTINInterpolator: public QgsInterpolator
{
public:
- QgsTINInterpolator( const QList<QgsVectorLayer*>& inputData );
+ QgsTINInterpolator( const QList<LayerData>& inputData, bool showProgressDialog = false );
~QgsTINInterpolator();
/**Calculates interpolation value for map coordinates x, y
@@ -37,12 +39,28 @@
@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
More information about the QGIS-commit
mailing list