From f351c4f7f8eaab4cc4a93a7a89a62783a45f0fe3 Mon Sep 17 00:00:00 2001 From: Marco Hugentobler Date: Thu, 11 Jul 2013 10:36:27 +0200 Subject: [PATCH] [FEATURE:] Add point sample class to analysis library --- python/analysis/analysis.sip | 1 + python/analysis/vector/qgspointsample.sip | 17 +++ src/analysis/CMakeLists.txt | 1 + src/analysis/vector/mersenne-twister.cpp | 1 - src/analysis/vector/mersenne-twister.h | 3 +- src/analysis/vector/qgspointsample.cpp | 163 ++++++++++++++++++++++ src/analysis/vector/qgspointsample.h | 42 ++++++ src/analysis/vector/qgstransectsample.cpp | 2 +- 8 files changed, 227 insertions(+), 3 deletions(-) create mode 100644 python/analysis/vector/qgspointsample.sip create mode 100644 src/analysis/vector/qgspointsample.cpp create mode 100644 src/analysis/vector/qgspointsample.h diff --git a/python/analysis/analysis.sip b/python/analysis/analysis.sip index e2409139174..3ced2905944 100644 --- a/python/analysis/analysis.sip +++ b/python/analysis/analysis.sip @@ -10,6 +10,7 @@ %Include vector/qgsgeometryanalyzer.sip %Include vector/qgsoverlayanalyzer.sip +%Include vector/qgspointsample.sip %Include vector/qgstransectsample.sip %Include vector/qgszonalstatistics.sip diff --git a/python/analysis/vector/qgspointsample.sip b/python/analysis/vector/qgspointsample.sip new file mode 100644 index 00000000000..8aa58cbe276 --- /dev/null +++ b/python/analysis/vector/qgspointsample.sip @@ -0,0 +1,17 @@ +/** \ingroup analysis + */ + +class QgsPointSample +{ +%TypeHeaderCode +#include +%End + + public: + QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute = QString() ); + ~QgsPointSample(); + + /**Starts calculation of random points + @return 0 in case of success*/ + int createRandomPoints( QProgressDialog* pd ); +}; diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index 8febaff63a5..b2bf6af8b35 100644 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -36,6 +36,7 @@ SET(QGIS_ANALYSIS_SRCS raster/qgsrastermatrix.cpp vector/mersenne-twister.cpp vector/qgsgeometryanalyzer.cpp + vector/qgspointsample.cpp vector/qgstransectsample.cpp vector/qgszonalstatistics.cpp vector/qgsoverlayanalyzer.cpp diff --git a/src/analysis/vector/mersenne-twister.cpp b/src/analysis/vector/mersenne-twister.cpp index e2d244d3287..dea4cb286d9 100644 --- a/src/analysis/vector/mersenne-twister.cpp +++ b/src/analysis/vector/mersenne-twister.cpp @@ -22,7 +22,6 @@ #include #include -#include #include "mersenne-twister.h" /* diff --git a/src/analysis/vector/mersenne-twister.h b/src/analysis/vector/mersenne-twister.h index c0ea845108c..49f76758a12 100644 --- a/src/analysis/vector/mersenne-twister.h +++ b/src/analysis/vector/mersenne-twister.h @@ -23,6 +23,7 @@ #define MERSENNE_TWISTER_H #include +#include #ifdef __cplusplus extern "C" { @@ -31,7 +32,7 @@ extern "C" { /* * Maximum number you can get from rand(). */ -#define RAND_MAX INT32_MAX +#define RAND_MAX std::numeric_limits::max() /* * Initialize the number generator with given seed. diff --git a/src/analysis/vector/qgspointsample.cpp b/src/analysis/vector/qgspointsample.cpp new file mode 100644 index 00000000000..82493ba3ebd --- /dev/null +++ b/src/analysis/vector/qgspointsample.cpp @@ -0,0 +1,163 @@ +#include "qgspointsample.h" +#include "qgsgeometry.h" +#include "qgsspatialindex.h" +#include "qgsvectorfilewriter.h" +#include "qgsvectorlayer.h" +#include +#include "mersenne-twister.h" + + +QgsPointSample::QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute ): mInputLayer( inputLayer ), + mOutputLayer( outputLayer ), mNumberOfPointsAttribute( nPointsAttribute ), mMinDistanceAttribute( minDistAttribute ), mNCreatedPoints( 0 ) +{ +} + +QgsPointSample::QgsPointSample() +{ +} + +QgsPointSample::~QgsPointSample() +{ +} + +int QgsPointSample::createRandomPoints( QProgressDialog* pd ) +{ + Q_UNUSED( pd ); + + //create input layer from id (test if polygon, valid) + if ( !mInputLayer ) + { + return 1; + } + + if ( mInputLayer->geometryType() != QGis::Polygon ) + { + return 2; + } + + //delete output file if it already exists + if ( QFile::exists( mOutputLayer ) ) + { + QgsVectorFileWriter::deleteShapeFile( mOutputLayer ); + } + + //create vector file writer + QgsFields outputFields; + outputFields.append( QgsField( "id", QVariant::Int ) ); + outputFields.append( QgsField( "station_id", QVariant::Int ) ); + outputFields.append( QgsField( "stratum_id", QVariant::Int ) ); + QgsVectorFileWriter writer( mOutputLayer, "UTF-8", + outputFields, + QGis::WKBPoint, + &( mInputLayer->crs() ) ); + + //check if creation of output layer successfull + if ( writer.hasError() != QgsVectorFileWriter::NoError ) + { + return 3; + } + + //init random number generator + mt_srand( QTime::currentTime().msec() ); + QgsFeature fet; + int nPoints = 0; + double minDistance = 0; + mNCreatedPoints = 0; + + QgsFeatureIterator fIt = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( + QStringList() << mNumberOfPointsAttribute << mMinDistanceAttribute, mInputLayer->pendingFields() ) ); + while ( fIt.nextFeature( fet ) ) + { + nPoints = fet.attribute( mNumberOfPointsAttribute ).toInt(); + if ( !mMinDistanceAttribute.isEmpty() ) + { + minDistance = fet.attribute( mMinDistanceAttribute ).toDouble(); + } + addSamplePoints( fet, writer, nPoints, minDistance ); + } + + return 0; +} + +void QgsPointSample::addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance ) +{ + QgsGeometry* geom = inputFeature.geometry(); + if ( !geom ) + { + return; + } + + QgsRectangle geomRect = geom->boundingBox(); + if ( geomRect.isEmpty() ) + { + return; + } + + QgsSpatialIndex sIndex; //to check minimum distance + QMap< QgsFeatureId, QgsPoint > pointMapForFeature; + + int nIterations = 0; + int maxIterations = nPoints * 200; + int points = 0; + + double randX = 0; + double randY = 0; + + while ( nIterations < maxIterations && points < nPoints ) + { + randX = (( double )mt_rand() / RAND_MAX ) * geomRect.width() + geomRect.xMinimum(); + randY = (( double )mt_rand() / RAND_MAX ) * geomRect.height() + geomRect.yMinimum(); + QgsPoint randPoint( randX, randY ); + QgsGeometry* ptGeom = QgsGeometry::fromPoint( randPoint ); + if ( ptGeom->within( geom ) && checkMinDistance( randPoint, sIndex, minDistance, pointMapForFeature ) ) + { + //add feature to writer + QgsFeature f( mNCreatedPoints ); + f.setAttribute( "id", mNCreatedPoints + 1 ); + f.setAttribute( "station_id", points + 1 ); + f.setAttribute( "stratum_id", inputFeature.id() ); + f.setGeometry( ptGeom ); + writer.addFeature( f ); + sIndex.insertFeature( f ); + pointMapForFeature.insert( mNCreatedPoints, randPoint ); + ++points; + ++mNCreatedPoints; + } + else + { + delete ptGeom; + } + ++nIterations; + } +} + +bool QgsPointSample::checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap ) +{ + if ( minDistance <= 0 ) + { + return true; + } + + QList neighborList = index.nearestNeighbor( pt, 1 ); + if ( neighborList.isEmpty() ) + { + return true; + } + + QMap< QgsFeatureId, QgsPoint >::const_iterator it = pointMap.find( neighborList[0] ); + if ( it == pointMap.constEnd() ) //should not happen + { + return true; + } + + QgsPoint neighborPt = it.value(); + if ( neighborPt.sqrDist( pt ) < ( minDistance * minDistance ) ) + { + return false; + } + return true; +} + + + + diff --git a/src/analysis/vector/qgspointsample.h b/src/analysis/vector/qgspointsample.h new file mode 100644 index 00000000000..7df28675faa --- /dev/null +++ b/src/analysis/vector/qgspointsample.h @@ -0,0 +1,42 @@ +#ifndef QGSPOINTSAMPLE_H +#define QGSPOINTSAMPLE_H + +#include "qgsfeature.h" +#include + +class QgsFeature; +class QgsPoint; +class QgsSpatialIndex; +class QgsVectorFileWriter; +class QgsVectorLayer; +class QProgressDialog; + +/**Creates random points in polygons / multipolygons*/ +class ANALYSIS_EXPORT QgsPointSample +{ + public: + QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute = QString() ); + ~QgsPointSample(); + + /**Starts calculation of random points + @return 0 in case of success*/ + int createRandomPoints( QProgressDialog* pd ); + + private: + + QgsPointSample(); //default constructor is forbidden + void addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance ); + bool checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap ); + + /**Layer id of input polygon/multipolygon layer*/ + QgsVectorLayer* mInputLayer; + /**Output path of result layer*/ + QString mOutputLayer; + /**Attribute containing number of points per feature*/ + QString mNumberOfPointsAttribute; + /**Attribute containing minimum distance between sample points (or -1 if no min. distance constraint)*/ + QString mMinDistanceAttribute; + QgsFeatureId mNCreatedPoints; //helper to find free ids +}; + +#endif // QGSPOINTSAMPLE_H diff --git a/src/analysis/vector/qgstransectsample.cpp b/src/analysis/vector/qgstransectsample.cpp index 2e6fb9316be..cd2d5ccf421 100644 --- a/src/analysis/vector/qgstransectsample.cpp +++ b/src/analysis/vector/qgstransectsample.cpp @@ -188,7 +188,7 @@ int QgsTransectSample::createSample( QProgressDialog* pd ) while ( nCreatedTransects < nTransects && nIterations < nMaxIterations ) { - double randomPosition = (( double )mt_rand() / std::numeric_limits::max() ) * clippedBaseline->length(); + double randomPosition = (( double )mt_rand() / RAND_MAX ) * clippedBaseline->length(); QgsGeometry* samplePoint = clippedBaseline->interpolate( randomPosition ); ++nIterations; if ( !samplePoint )