Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmdbscanclustering.cpp
3 : : ---------------------
4 : : begin : July 2018
5 : : copyright : (C) 2018 by Nyall Dawson
6 : : email : nyall dot dawson at gmail dot com
7 : : ***************************************************************************/
8 : :
9 : : /***************************************************************************
10 : : * *
11 : : * This program is free software; you can redistribute it and/or modify *
12 : : * it under the terms of the GNU General Public License as published by *
13 : : * the Free Software Foundation; either version 2 of the License, or *
14 : : * (at your option) any later version. *
15 : : * *
16 : : ***************************************************************************/
17 : :
18 : : #include "qgsalgorithmdbscanclustering.h"
19 : : #include "qgsspatialindexkdbush.h"
20 : : #include <unordered_set>
21 : :
22 : : ///@cond PRIVATE
23 : :
24 : 0 : QString QgsDbscanClusteringAlgorithm::name() const
25 : : {
26 : 0 : return QStringLiteral( "dbscanclustering" );
27 : : }
28 : :
29 : 0 : QString QgsDbscanClusteringAlgorithm::displayName() const
30 : : {
31 : 0 : return QObject::tr( "DBSCAN clustering" );
32 : : }
33 : :
34 : 0 : QString QgsDbscanClusteringAlgorithm::shortDescription() const
35 : : {
36 : 0 : return QObject::tr( "Clusters point features using a density based scan algorithm." );
37 : : }
38 : :
39 : 0 : QStringList QgsDbscanClusteringAlgorithm::tags() const
40 : : {
41 : 0 : return QObject::tr( "clustering,clusters,density,based,points" ).split( ',' );
42 : 0 : }
43 : :
44 : 0 : QString QgsDbscanClusteringAlgorithm::group() const
45 : : {
46 : 0 : return QObject::tr( "Vector analysis" );
47 : : }
48 : :
49 : 0 : QString QgsDbscanClusteringAlgorithm::groupId() const
50 : : {
51 : 0 : return QStringLiteral( "vectoranalysis" );
52 : : }
53 : :
54 : 0 : void QgsDbscanClusteringAlgorithm::initAlgorithm( const QVariantMap & )
55 : : {
56 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ),
57 : 0 : QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorPoint ) );
58 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MIN_SIZE" ), QObject::tr( "Minimum cluster size" ),
59 : 0 : QgsProcessingParameterNumber::Integer, 5, false, 1 ) );
60 : 0 : addParameter( new QgsProcessingParameterDistance( QStringLiteral( "EPS" ),
61 : 0 : QObject::tr( "Maximum distance between clustered points" ), 1, QStringLiteral( "INPUT" ), false, 0 ) );
62 : :
63 : 0 : auto dbscanStarParam = std::make_unique<QgsProcessingParameterBoolean>( QStringLiteral( "DBSCAN*" ),
64 : 0 : QObject::tr( "Treat border points as noise (DBSCAN*)" ), false, true );
65 : 0 : dbscanStarParam->setFlags( dbscanStarParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
66 : 0 : addParameter( dbscanStarParam.release() );
67 : :
68 : 0 : auto fieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "FIELD_NAME" ),
69 : 0 : QObject::tr( "Cluster field name" ), QStringLiteral( "CLUSTER_ID" ) );
70 : 0 : fieldNameParam->setFlags( fieldNameParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
71 : 0 : addParameter( fieldNameParam.release() );
72 : 0 : auto sizeFieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "SIZE_FIELD_NAME" ),
73 : 0 : QObject::tr( "Cluster size field name" ), QStringLiteral( "CLUSTER_SIZE" ) );
74 : 0 : sizeFieldNameParam->setFlags( sizeFieldNameParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
75 : 0 : addParameter( sizeFieldNameParam.release() );
76 : :
77 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Clusters" ), QgsProcessing::TypeVectorPoint ) );
78 : :
79 : 0 : addOutput( new QgsProcessingOutputNumber( QStringLiteral( "NUM_CLUSTERS" ), QObject::tr( "Number of clusters" ) ) );
80 : 0 : }
81 : :
82 : 0 : QString QgsDbscanClusteringAlgorithm::shortHelpString() const
83 : : {
84 : 0 : return QObject::tr( "Clusters point features based on a 2D implementation of Density-based spatial clustering of applications with noise (DBSCAN) algorithm.\n\n"
85 : : "The algorithm requires two parameters, a minimum cluster size (“minPts”), and the maximum distance allowed between clustered points (“eps”)." );
86 : : }
87 : :
88 : 0 : QgsDbscanClusteringAlgorithm *QgsDbscanClusteringAlgorithm::createInstance() const
89 : : {
90 : 0 : return new QgsDbscanClusteringAlgorithm();
91 : : }
92 : :
93 : : struct KDBushDataEqualById
94 : : {
95 : 0 : bool operator()( const QgsSpatialIndexKDBushData &a, const QgsSpatialIndexKDBushData &b ) const
96 : : {
97 : 0 : return a.id == b.id;
98 : : }
99 : : };
100 : :
101 : : struct KDBushDataHashById
102 : : {
103 : 0 : std::size_t operator()( const QgsSpatialIndexKDBushData &a ) const
104 : : {
105 : 0 : return std::hash< QgsFeatureId > {}( a.id );
106 : : }
107 : : };
108 : :
109 : 0 : QVariantMap QgsDbscanClusteringAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
110 : : {
111 : 0 : std::unique_ptr< QgsProcessingFeatureSource > source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
112 : 0 : if ( !source )
113 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
114 : :
115 : 0 : const std::size_t minSize = static_cast< std::size_t>( parameterAsInt( parameters, QStringLiteral( "MIN_SIZE" ), context ) );
116 : 0 : const double eps = parameterAsDouble( parameters, QStringLiteral( "EPS" ), context );
117 : 0 : const bool borderPointsAreNoise = parameterAsBoolean( parameters, QStringLiteral( "DBSCAN*" ), context );
118 : :
119 : 0 : QgsFields outputFields = source->fields();
120 : 0 : QgsFields newFields;
121 : 0 : const QString clusterFieldName = parameterAsString( parameters, QStringLiteral( "FIELD_NAME" ), context );
122 : 0 : newFields.append( QgsField( clusterFieldName, QVariant::Int ) );
123 : 0 : const QString clusterSizeFieldName = parameterAsString( parameters, QStringLiteral( "SIZE_FIELD_NAME" ), context );
124 : 0 : newFields.append( QgsField( clusterSizeFieldName, QVariant::Int ) );
125 : 0 : outputFields = QgsProcessingUtils::combineFields( outputFields, newFields );
126 : :
127 : 0 : QString dest;
128 : 0 : std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, outputFields, source->wkbType(), source->sourceCrs() ) );
129 : 0 : if ( !sink )
130 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
131 : :
132 : : // build spatial index
133 : 0 : feedback->pushInfo( QObject::tr( "Building spatial index" ) );
134 : 0 : QgsSpatialIndexKDBush index( *source, feedback );
135 : 0 : if ( feedback->isCanceled() )
136 : 0 : return QVariantMap();
137 : :
138 : : // dbscan!
139 : 0 : feedback->pushInfo( QObject::tr( "Analysing clusters" ) );
140 : 0 : std::unordered_map< QgsFeatureId, int> idToCluster;
141 : 0 : idToCluster.reserve( index.size() );
142 : 0 : QgsFeatureIterator features = source->getFeatures( QgsFeatureRequest().setNoAttributes() );
143 : 0 : const long featureCount = source->featureCount();
144 : 0 : dbscan( minSize, eps, borderPointsAreNoise, featureCount, features, index, idToCluster, feedback );
145 : :
146 : : // cluster size
147 : 0 : std::unordered_map< int, int> clusterSize;
148 : 0 : std::for_each( idToCluster.begin(), idToCluster.end(), [ &clusterSize ]( std::pair< QgsFeatureId, int > idCluster ) { clusterSize[ idCluster.second ]++; } );
149 : :
150 : : // write clusters
151 : 0 : const double writeStep = featureCount > 0 ? 10.0 / featureCount : 1;
152 : 0 : features = source->getFeatures();
153 : 0 : int i = 0;
154 : 0 : QgsFeature feat;
155 : 0 : while ( features.nextFeature( feat ) )
156 : : {
157 : 0 : i++;
158 : 0 : if ( feedback->isCanceled() )
159 : : {
160 : 0 : break;
161 : : }
162 : :
163 : 0 : feedback->setProgress( 90 + i * writeStep );
164 : 0 : QgsAttributes attr = feat.attributes();
165 : 0 : auto cluster = idToCluster.find( feat.id() );
166 : 0 : if ( cluster != idToCluster.end() )
167 : : {
168 : 0 : attr << cluster->second << clusterSize[ cluster->second ];
169 : 0 : }
170 : : else
171 : : {
172 : 0 : attr << QVariant() << QVariant();
173 : : }
174 : 0 : feat.setAttributes( attr );
175 : 0 : sink->addFeature( feat, QgsFeatureSink::FastInsert );
176 : 0 : }
177 : :
178 : 0 : QVariantMap outputs;
179 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), dest );
180 : 0 : outputs.insert( QStringLiteral( "NUM_CLUSTERS" ), static_cast< unsigned int >( clusterSize.size() ) );
181 : 0 : return outputs;
182 : 0 : }
183 : :
184 : :
185 : 0 : void QgsDbscanClusteringAlgorithm::dbscan( const std::size_t minSize,
186 : : const double eps,
187 : : const bool borderPointsAreNoise,
188 : : const long featureCount,
189 : : QgsFeatureIterator features,
190 : : QgsSpatialIndexKDBush &index,
191 : : std::unordered_map< QgsFeatureId, int> &idToCluster,
192 : : QgsProcessingFeedback *feedback )
193 : : {
194 : 0 : const double step = featureCount > 0 ? 90.0 / featureCount : 1;
195 : :
196 : 0 : std::unordered_set< QgsFeatureId > visited;
197 : 0 : visited.reserve( index.size() );
198 : :
199 : 0 : QgsFeature feat;
200 : 0 : int i = 0;
201 : 0 : int clusterCount = 0;
202 : :
203 : 0 : while ( features.nextFeature( feat ) )
204 : : {
205 : 0 : if ( feedback->isCanceled() )
206 : : {
207 : 0 : break;
208 : : }
209 : :
210 : 0 : if ( !feat.hasGeometry() )
211 : : {
212 : 0 : feedback->setProgress( ++i * step );
213 : 0 : continue;
214 : : }
215 : :
216 : 0 : if ( visited.find( feat.id() ) != visited.end() )
217 : : {
218 : : // already visited!
219 : 0 : continue;
220 : : }
221 : :
222 : 0 : QgsPointXY point;
223 : 0 : if ( QgsWkbTypes::flatType( feat.geometry().wkbType() ) == QgsWkbTypes::Point )
224 : 0 : point = QgsPointXY( *qgsgeometry_cast< const QgsPoint * >( feat.geometry().constGet() ) );
225 : : else
226 : : {
227 : : // not a point geometry
228 : 0 : feedback->reportError( QObject::tr( "Feature %1 is a %2 feature, not a point." ).arg( feat.id() ).arg( QgsWkbTypes::displayString( feat.geometry().wkbType() ) ) );
229 : 0 : feedback->setProgress( ++i * step );
230 : 0 : continue;
231 : : }
232 : :
233 : 0 : std::unordered_set< QgsSpatialIndexKDBushData, KDBushDataHashById, KDBushDataEqualById> within;
234 : :
235 : 0 : if ( minSize > 1 )
236 : : {
237 : 0 : index.within( point, eps, [ &within]( const QgsSpatialIndexKDBushData & data )
238 : : {
239 : 0 : within.insert( data );
240 : 0 : } );
241 : 0 : if ( within.size() < minSize )
242 : 0 : continue;
243 : :
244 : 0 : visited.insert( feat.id() );
245 : 0 : }
246 : : else
247 : : {
248 : : // optimised case for minSize == 1, we can skip the initial check
249 : 0 : within.insert( QgsSpatialIndexKDBushData( feat.id(), point.x(), point.y() ) );
250 : : }
251 : :
252 : : // start new cluster
253 : 0 : clusterCount++;
254 : 0 : idToCluster[ feat.id() ] = clusterCount;
255 : 0 : feedback->setProgress( ++i * step );
256 : :
257 : 0 : while ( !within.empty() )
258 : : {
259 : 0 : if ( feedback->isCanceled() )
260 : : {
261 : 0 : break;
262 : : }
263 : :
264 : 0 : QgsSpatialIndexKDBushData j = *within.begin();
265 : 0 : within.erase( within.begin() );
266 : :
267 : 0 : if ( visited.find( j.id ) != visited.end() )
268 : : {
269 : : // already visited!
270 : 0 : continue;
271 : : }
272 : :
273 : 0 : visited.insert( j.id );
274 : 0 : feedback->setProgress( ++i * step );
275 : :
276 : : // check from this point
277 : 0 : QgsPointXY point2 = j.point();
278 : :
279 : 0 : std::unordered_set< QgsSpatialIndexKDBushData, KDBushDataHashById, KDBushDataEqualById > within2;
280 : 0 : index.within( point2, eps, [&within2]( const QgsSpatialIndexKDBushData & data )
281 : : {
282 : 0 : within2.insert( data );
283 : 0 : } );
284 : 0 : if ( within2.size() >= minSize )
285 : : {
286 : : // expand neighbourhood
287 : 0 : std::copy_if( within2.begin(),
288 : 0 : within2.end(),
289 : 0 : std::inserter( within, within.end() ),
290 : 0 : [&visited]( const QgsSpatialIndexKDBushData & needle )
291 : : {
292 : 0 : return visited.find( needle.id ) == visited.end();
293 : : } );
294 : 0 : }
295 : 0 : if ( !borderPointsAreNoise || within2.size() >= minSize )
296 : : {
297 : 0 : idToCluster[ j.id ] = clusterCount;
298 : 0 : }
299 : 0 : }
300 : 0 : }
301 : 0 : }
302 : :
303 : : ///@endcond
304 : :
305 : :
|
