QGIS API Documentation 3.41.0-Master (d5b93354e9c)
Loading...
Searching...
No Matches
qgsalgorithmkmeansclustering.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmkmeansclustering.cpp
3 ---------------------
4 begin : June 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
19#include <unordered_map>
20
22
23const int KMEANS_MAX_ITERATIONS = 1000;
24
25QString QgsKMeansClusteringAlgorithm::name() const
26{
27 return QStringLiteral( "kmeansclustering" );
28}
29
30QString QgsKMeansClusteringAlgorithm::displayName() const
31{
32 return QObject::tr( "K-means clustering" );
33}
34
35QStringList QgsKMeansClusteringAlgorithm::tags() const
36{
37 return QObject::tr( "clustering,clusters,kmeans,points" ).split( ',' );
38}
39
40QString QgsKMeansClusteringAlgorithm::group() const
41{
42 return QObject::tr( "Vector analysis" );
43}
44
45QString QgsKMeansClusteringAlgorithm::groupId() const
46{
47 return QStringLiteral( "vectoranalysis" );
48}
49
50void QgsKMeansClusteringAlgorithm::initAlgorithm( const QVariantMap & )
51{
52 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorAnyGeometry ) ) );
53 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "CLUSTERS" ), QObject::tr( "Number of clusters" ), Qgis::ProcessingNumberParameterType::Integer, 5, false, 1 ) );
54
55 auto fieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "FIELD_NAME" ), QObject::tr( "Cluster field name" ), QStringLiteral( "CLUSTER_ID" ) );
56 fieldNameParam->setFlags( fieldNameParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
57 addParameter( fieldNameParam.release() );
58 auto sizeFieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "SIZE_FIELD_NAME" ), QObject::tr( "Cluster size field name" ), QStringLiteral( "CLUSTER_SIZE" ) );
59 sizeFieldNameParam->setFlags( sizeFieldNameParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
60 addParameter( sizeFieldNameParam.release() );
61
62 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Clusters" ), Qgis::ProcessingSourceType::VectorAnyGeometry ) );
63}
64
65QString QgsKMeansClusteringAlgorithm::shortHelpString() const
66{
67 return QObject::tr( "Calculates the 2D distance based k-means cluster number for each input feature.\n\n"
68 "If input geometries are lines or polygons, the clustering is based on the centroid of the feature." );
69}
70
71QgsKMeansClusteringAlgorithm *QgsKMeansClusteringAlgorithm::createInstance() const
72{
73 return new QgsKMeansClusteringAlgorithm();
74}
75
76QVariantMap QgsKMeansClusteringAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
77{
78 std::unique_ptr<QgsProcessingFeatureSource> source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
79 if ( !source )
80 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
81
82 int k = parameterAsInt( parameters, QStringLiteral( "CLUSTERS" ), context );
83
84 QgsFields outputFields = source->fields();
85 QgsFields newFields;
86 const QString clusterFieldName = parameterAsString( parameters, QStringLiteral( "FIELD_NAME" ), context );
87 newFields.append( QgsField( clusterFieldName, QMetaType::Type::Int ) );
88 const QString clusterSizeFieldName = parameterAsString( parameters, QStringLiteral( "SIZE_FIELD_NAME" ), context );
89 newFields.append( QgsField( clusterSizeFieldName, QMetaType::Type::Int ) );
90 outputFields = QgsProcessingUtils::combineFields( outputFields, newFields );
91
92 QString dest;
93 std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, outputFields, source->wkbType(), source->sourceCrs() ) );
94 if ( !sink )
95 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
96
97 // build list of point inputs - if it's already a point, use that. If not, take the centroid.
98 feedback->pushInfo( QObject::tr( "Collecting input points" ) );
99 const double step = source->featureCount() > 0 ? 50.0 / source->featureCount() : 1;
100 int i = 0;
101 int n = 0;
102 int featureWithGeometryCount = 0;
103 QgsFeature feat;
104
105 std::vector<Feature> clusterFeatures;
106 QgsFeatureIterator features = source->getFeatures( QgsFeatureRequest().setNoAttributes() );
107 QHash<QgsFeatureId, int> idToObj;
108 while ( features.nextFeature( feat ) )
109 {
110 i++;
111 if ( feedback->isCanceled() )
112 {
113 break;
114 }
115
116 feedback->setProgress( i * step );
117 if ( !feat.hasGeometry() )
118 continue;
119 featureWithGeometryCount++;
120
121 QgsPointXY point;
123 point = QgsPointXY( *qgsgeometry_cast<const QgsPoint *>( feat.geometry().constGet() ) );
124 else
125 {
126 const QgsGeometry centroid = feat.geometry().centroid();
127 if ( centroid.isNull() )
128 continue; // centroid failed, e.g. empty linestring
129
130 point = QgsPointXY( *qgsgeometry_cast<const QgsPoint *>( centroid.constGet() ) );
131 }
132
133 n++;
134
135 idToObj[feat.id()] = clusterFeatures.size();
136 clusterFeatures.emplace_back( Feature( point ) );
137 }
138
139 if ( n < k )
140 {
141 feedback->reportError( QObject::tr( "Number of geometries is less than the number of clusters requested, not all clusters will get data" ) );
142 k = n;
143 }
144
145 if ( k > 1 )
146 {
147 feedback->pushInfo( QObject::tr( "Calculating clusters" ) );
148
149 // cluster centers
150 std::vector<QgsPointXY> centers( k );
151
152 initClusters( clusterFeatures, centers, k, feedback );
153 calculateKMeans( clusterFeatures, centers, k, feedback );
154 }
155
156 // cluster size
157 std::unordered_map<int, int> clusterSize;
158 for ( const int obj : idToObj )
159 clusterSize[clusterFeatures[obj].cluster]++;
160
161 features = source->getFeatures();
162 i = 0;
163 while ( features.nextFeature( feat ) )
164 {
165 i++;
166 if ( feedback->isCanceled() )
167 {
168 break;
169 }
170
171 feedback->setProgress( 50 + i * step );
172 QgsAttributes attr = feat.attributes();
173 const auto obj = idToObj.find( feat.id() );
174 if ( !feat.hasGeometry() || obj == idToObj.end() )
175 {
176 attr << QVariant() << QVariant();
177 }
178 else if ( k <= 1 )
179 {
180 attr << 0 << featureWithGeometryCount;
181 }
182 else
183 {
184 const int cluster = clusterFeatures[*obj].cluster;
185 attr << cluster << clusterSize[cluster];
186 }
187 feat.setAttributes( attr );
188 if ( !sink->addFeature( feat, QgsFeatureSink::FastInsert ) )
189 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
190 }
191
192 sink->finalize();
193
194 QVariantMap outputs;
195 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
196 return outputs;
197}
198
199// ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
200
201void QgsKMeansClusteringAlgorithm::initClusters( std::vector<Feature> &points, std::vector<QgsPointXY> &centers, const int k, QgsProcessingFeedback *feedback )
202{
203 const std::size_t n = points.size();
204 if ( n == 0 )
205 return;
206
207 if ( n == 1 )
208 {
209 for ( int i = 0; i < k; i++ )
210 centers[i] = points[0].point;
211 return;
212 }
213
214 long duplicateCount = 1;
215 // initially scan for two most distance points from each other, p1 and p2
216 std::size_t p1 = 0;
217 std::size_t p2 = 0;
218 double distanceP1 = 0;
219 double distanceP2 = 0;
220 double maxDistance = -1;
221 for ( std::size_t i = 1; i < n; i++ )
222 {
223 distanceP1 = points[i].point.sqrDist( points[p1].point );
224 distanceP2 = points[i].point.sqrDist( points[p2].point );
225
226 // if this point is further then existing candidates, replace our choice
227 if ( ( distanceP1 > maxDistance ) || ( distanceP2 > maxDistance ) )
228 {
229 maxDistance = std::max( distanceP1, distanceP2 );
230 if ( distanceP1 > distanceP2 )
231 p2 = i;
232 else
233 p1 = i;
234 }
235
236 // also record count of duplicate points
237 if ( qgsDoubleNear( distanceP1, 0 ) || qgsDoubleNear( distanceP2, 0 ) )
238 duplicateCount++;
239 }
240
241 if ( feedback && duplicateCount > 1 )
242 {
243 feedback->pushInfo( QObject::tr( "There are at least %n duplicate input(s), the number of output clusters may be less than was requested", nullptr, duplicateCount ) );
244 }
245
246 // By now two points should be found and be not the same
247 // Q_ASSERT( p1 != p2 && maxDistance >= 0 );
248
249 // Accept these two points as our initial centers
250 centers[0] = points[p1].point;
251 centers[1] = points[p2].point;
252
253 if ( k > 2 )
254 {
255 // array of minimum distance to a point from accepted cluster centers
256 std::vector<double> distances( n );
257
258 // initialize array with distance to first object
259 for ( std::size_t j = 0; j < n; j++ )
260 {
261 distances[j] = points[j].point.sqrDist( centers[0] );
262 }
263 distances[p1] = -1;
264 distances[p2] = -1;
265
266 // loop i on clusters, skip 0 and 1 as found already
267 for ( int i = 2; i < k; i++ )
268 {
269 std::size_t candidateCenter = 0;
270 double maxDistance = std::numeric_limits<double>::lowest();
271
272 // loop j on points
273 for ( std::size_t j = 0; j < n; j++ )
274 {
275 // accepted clusters are already marked with distance = -1
276 if ( distances[j] < 0 )
277 continue;
278
279 // update minimal distance with previously accepted cluster
280 distances[j] = std::min( points[j].point.sqrDist( centers[i - 1] ), distances[j] );
281
282 // greedily take a point that's farthest from any of accepted clusters
283 if ( distances[j] > maxDistance )
284 {
285 candidateCenter = j;
286 maxDistance = distances[j];
287 }
288 }
289
290 // checked earlier by counting entries on input, just in case
291 Q_ASSERT( maxDistance >= 0 );
292
293 // accept candidate to centers
294 distances[candidateCenter] = -1;
295 // copy the point coordinates into the initial centers array
296 centers[i] = points[candidateCenter].point;
297 }
298 }
299}
300
301// ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
302
303void QgsKMeansClusteringAlgorithm::calculateKMeans( std::vector<QgsKMeansClusteringAlgorithm::Feature> &objs, std::vector<QgsPointXY> &centers, int k, QgsProcessingFeedback *feedback )
304{
305 int converged = false;
306 bool changed = false;
307
308 // avoid reallocating weights array for every iteration
309 std::vector<uint> weights( k );
310
311 uint i = 0;
312 for ( i = 0; i < KMEANS_MAX_ITERATIONS && !converged; i++ )
313 {
314 if ( feedback && feedback->isCanceled() )
315 break;
316
317 findNearest( objs, centers, k, changed );
318 updateMeans( objs, centers, weights, k );
319 converged = !changed;
320 }
321
322 if ( !converged && feedback )
323 feedback->reportError( QObject::tr( "Clustering did not converge after %n iteration(s)", nullptr, i ) );
324 else if ( feedback )
325 feedback->pushInfo( QObject::tr( "Clustering converged after %n iteration(s)", nullptr, i ) );
326}
327
328// ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
329
330void QgsKMeansClusteringAlgorithm::findNearest( std::vector<QgsKMeansClusteringAlgorithm::Feature> &points, const std::vector<QgsPointXY> &centers, const int k, bool &changed )
331{
332 changed = false;
333 const std::size_t n = points.size();
334 for ( std::size_t i = 0; i < n; i++ )
335 {
336 Feature &point = points[i];
337
338 // Initialize with distance to first cluster
339 double currentDistance = point.point.sqrDist( centers[0] );
340 int currentCluster = 0;
341
342 // Check all other cluster centers and find the nearest
343 for ( int cluster = 1; cluster < k; cluster++ )
344 {
345 const double distance = point.point.sqrDist( centers[cluster] );
346 if ( distance < currentDistance )
347 {
348 currentDistance = distance;
349 currentCluster = cluster;
350 }
351 }
352
353 // Store the nearest cluster this object is in
354 if ( point.cluster != currentCluster )
355 {
356 changed = true;
357 point.cluster = currentCluster;
358 }
359 }
360}
361
362// ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
363
364void QgsKMeansClusteringAlgorithm::updateMeans( const std::vector<Feature> &points, std::vector<QgsPointXY> &centers, std::vector<uint> &weights, const int k )
365{
366 const uint n = points.size();
367 std::fill( weights.begin(), weights.end(), 0 );
368 for ( int i = 0; i < k; i++ )
369 {
370 centers[i].setX( 0.0 );
371 centers[i].setY( 0.0 );
372 }
373 for ( uint i = 0; i < n; i++ )
374 {
375 const int cluster = points[i].cluster;
376 centers[cluster] += QgsVector( points[i].point.x(), points[i].point.y() );
377 weights[cluster] += 1;
378 }
379 for ( int i = 0; i < k; i++ )
380 {
381 centers[i] /= weights[i];
382 }
383}
384
385
@ VectorAnyGeometry
Any vector layer with geometry.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
A vector of attributes.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsAttributes attributes
Definition qgsfeature.h:67
QgsFeatureId id
Definition qgsfeature.h:66
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:70
A geometry is the spatial representation of a feature.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsGeometry centroid() const
Returns the center of mass of a geometry.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
A class to represent a 2D point.
Definition qgspointxy.h:60
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
virtual void pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A numeric parameter for processing algorithms.
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
A class to represent a vector.
Definition qgsvector.h:30
static Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB type.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:6024