Branch data Line data Source code
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 : :
18 : : #include "qgsalgorithmkmeansclustering.h"
19 : : #include <unordered_map>
20 : :
21 : : ///@cond PRIVATE
22 : :
23 : : const int KMEANS_MAX_ITERATIONS = 1000;
24 : :
25 : 0 : QString QgsKMeansClusteringAlgorithm::name() const
26 : : {
27 : 0 : return QStringLiteral( "kmeansclustering" );
28 : : }
29 : :
30 : 0 : QString QgsKMeansClusteringAlgorithm::displayName() const
31 : : {
32 : 0 : return QObject::tr( "K-means clustering" );
33 : : }
34 : :
35 : 0 : QStringList QgsKMeansClusteringAlgorithm::tags() const
36 : : {
37 : 0 : return QObject::tr( "clustering,clusters,kmeans,points" ).split( ',' );
38 : 0 : }
39 : :
40 : 0 : QString QgsKMeansClusteringAlgorithm::group() const
41 : : {
42 : 0 : return QObject::tr( "Vector analysis" );
43 : : }
44 : :
45 : 0 : QString QgsKMeansClusteringAlgorithm::groupId() const
46 : : {
47 : 0 : return QStringLiteral( "vectoranalysis" );
48 : : }
49 : :
50 : 0 : void QgsKMeansClusteringAlgorithm::initAlgorithm( const QVariantMap & )
51 : : {
52 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ),
53 : 0 : QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorAnyGeometry ) );
54 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "CLUSTERS" ), QObject::tr( "Number of clusters" ),
55 : 0 : QgsProcessingParameterNumber::Integer, 5, false, 1 ) );
56 : :
57 : 0 : auto fieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "FIELD_NAME" ),
58 : 0 : QObject::tr( "Cluster field name" ), QStringLiteral( "CLUSTER_ID" ) );
59 : 0 : fieldNameParam->setFlags( fieldNameParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
60 : 0 : addParameter( fieldNameParam.release() );
61 : 0 : auto sizeFieldNameParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "SIZE_FIELD_NAME" ),
62 : 0 : QObject::tr( "Cluster size field name" ), QStringLiteral( "CLUSTER_SIZE" ) );
63 : 0 : sizeFieldNameParam->setFlags( sizeFieldNameParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
64 : 0 : addParameter( sizeFieldNameParam.release() );
65 : :
66 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Clusters" ), QgsProcessing::TypeVectorAnyGeometry ) );
67 : 0 : }
68 : :
69 : 0 : QString QgsKMeansClusteringAlgorithm::shortHelpString() const
70 : : {
71 : 0 : return QObject::tr( "Calculates the 2D distance based k-means cluster number for each input feature.\n\n"
72 : : "If input geometries are lines or polygons, the clustering is based on the centroid of the feature." );
73 : : }
74 : :
75 : 0 : QgsKMeansClusteringAlgorithm *QgsKMeansClusteringAlgorithm::createInstance() const
76 : : {
77 : 0 : return new QgsKMeansClusteringAlgorithm();
78 : : }
79 : :
80 : 0 : QVariantMap QgsKMeansClusteringAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
81 : : {
82 : 0 : std::unique_ptr< QgsProcessingFeatureSource > source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
83 : 0 : if ( !source )
84 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
85 : :
86 : 0 : int k = parameterAsInt( parameters, QStringLiteral( "CLUSTERS" ), context );
87 : :
88 : 0 : QgsFields outputFields = source->fields();
89 : 0 : QgsFields newFields;
90 : 0 : const QString clusterFieldName = parameterAsString( parameters, QStringLiteral( "FIELD_NAME" ), context );
91 : 0 : newFields.append( QgsField( clusterFieldName, QVariant::Int ) );
92 : 0 : const QString clusterSizeFieldName = parameterAsString( parameters, QStringLiteral( "SIZE_FIELD_NAME" ), context );
93 : 0 : newFields.append( QgsField( clusterSizeFieldName, QVariant::Int ) );
94 : 0 : outputFields = QgsProcessingUtils::combineFields( outputFields, newFields );
95 : :
96 : 0 : QString dest;
97 : 0 : std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, outputFields, source->wkbType(), source->sourceCrs() ) );
98 : 0 : if ( !sink )
99 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
100 : :
101 : : // build list of point inputs - if it's already a point, use that. If not, take the centroid.
102 : 0 : feedback->pushInfo( QObject::tr( "Collecting input points" ) );
103 : 0 : double step = source->featureCount() > 0 ? 50.0 / source->featureCount() : 1;
104 : 0 : int i = 0;
105 : 0 : int n = 0;
106 : 0 : int featureWithGeometryCount = 0;
107 : 0 : QgsFeature feat;
108 : :
109 : 0 : std::vector< Feature > clusterFeatures;
110 : 0 : QgsFeatureIterator features = source->getFeatures( QgsFeatureRequest().setNoAttributes() );
111 : 0 : QHash< QgsFeatureId, int > idToObj;
112 : 0 : while ( features.nextFeature( feat ) )
113 : : {
114 : 0 : i++;
115 : 0 : if ( feedback->isCanceled() )
116 : : {
117 : 0 : break;
118 : : }
119 : :
120 : 0 : feedback->setProgress( i * step );
121 : 0 : if ( !feat.hasGeometry() )
122 : 0 : continue;
123 : 0 : featureWithGeometryCount++;
124 : :
125 : 0 : QgsPointXY point;
126 : 0 : if ( QgsWkbTypes::flatType( feat.geometry().wkbType() ) == QgsWkbTypes::Point )
127 : 0 : point = QgsPointXY( *qgsgeometry_cast< const QgsPoint * >( feat.geometry().constGet() ) );
128 : : else
129 : : {
130 : 0 : QgsGeometry centroid = feat.geometry().centroid();
131 : 0 : if ( centroid.isNull() )
132 : 0 : continue; // centroid failed, e.g. empty linestring
133 : :
134 : 0 : point = QgsPointXY( *qgsgeometry_cast< const QgsPoint * >( centroid.constGet() ) );
135 : 0 : }
136 : :
137 : 0 : n++;
138 : :
139 : 0 : idToObj[ feat.id() ] = clusterFeatures.size();
140 : 0 : clusterFeatures.emplace_back( Feature( point ) );
141 : : }
142 : :
143 : 0 : if ( n < k )
144 : : {
145 : 0 : feedback->reportError( QObject::tr( "Number of geometries is less than the number of clusters requested, not all clusters will get data" ) );
146 : 0 : k = n;
147 : 0 : }
148 : :
149 : 0 : if ( k > 1 )
150 : : {
151 : 0 : feedback->pushInfo( QObject::tr( "Calculating clusters" ) );
152 : :
153 : : // cluster centers
154 : 0 : std::vector< QgsPointXY > centers( k );
155 : :
156 : 0 : initClusters( clusterFeatures, centers, k, feedback );
157 : 0 : calculateKMeans( clusterFeatures, centers, k, feedback );
158 : 0 : }
159 : :
160 : : // cluster size
161 : 0 : std::unordered_map< int, int> clusterSize;
162 : 0 : for ( int obj : idToObj )
163 : 0 : clusterSize[ clusterFeatures[ obj ].cluster ]++;
164 : :
165 : 0 : features = source->getFeatures();
166 : 0 : i = 0;
167 : 0 : while ( features.nextFeature( feat ) )
168 : : {
169 : 0 : i++;
170 : 0 : if ( feedback->isCanceled() )
171 : : {
172 : 0 : break;
173 : : }
174 : :
175 : 0 : feedback->setProgress( 50 + i * step );
176 : 0 : QgsAttributes attr = feat.attributes();
177 : 0 : auto obj = idToObj.find( feat.id() );
178 : 0 : if ( !feat.hasGeometry() || obj == idToObj.end() )
179 : : {
180 : 0 : attr << QVariant() << QVariant();
181 : 0 : }
182 : 0 : else if ( k <= 1 )
183 : : {
184 : 0 : attr << 0 << featureWithGeometryCount;
185 : 0 : }
186 : : else
187 : : {
188 : 0 : int cluster = clusterFeatures[ *obj ].cluster;
189 : 0 : attr << cluster << clusterSize[ cluster ];
190 : : }
191 : 0 : feat.setAttributes( attr );
192 : 0 : sink->addFeature( feat, QgsFeatureSink::FastInsert );
193 : 0 : }
194 : :
195 : 0 : QVariantMap outputs;
196 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), dest );
197 : 0 : return outputs;
198 : 0 : }
199 : :
200 : : // ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
201 : :
202 : 0 : void QgsKMeansClusteringAlgorithm::initClusters( std::vector<Feature> &points, std::vector<QgsPointXY> ¢ers, const int k, QgsProcessingFeedback *feedback )
203 : : {
204 : 0 : std::size_t n = points.size();
205 : 0 : if ( n == 0 )
206 : 0 : return;
207 : :
208 : 0 : if ( n == 1 )
209 : : {
210 : 0 : for ( int i = 0; i < k; i++ )
211 : 0 : centers[ i ] = points[ 0 ].point;
212 : 0 : return;
213 : : }
214 : :
215 : 0 : long duplicateCount = 1;
216 : : // initially scan for two most distance points from each other, p1 and p2
217 : 0 : std::size_t p1 = 0;
218 : 0 : std::size_t p2 = 0;
219 : 0 : double distanceP1 = 0;
220 : 0 : double distanceP2 = 0;
221 : 0 : double maxDistance = -1;
222 : 0 : for ( std::size_t i = 1; i < n; i++ )
223 : : {
224 : 0 : distanceP1 = points[i].point.sqrDist( points[p1].point );
225 : 0 : distanceP2 = points[i].point.sqrDist( points[p2].point );
226 : :
227 : : // if this point is further then existing candidates, replace our choice
228 : 0 : if ( ( distanceP1 > maxDistance ) || ( distanceP2 > maxDistance ) )
229 : : {
230 : 0 : maxDistance = std::max( distanceP1, distanceP2 );
231 : 0 : if ( distanceP1 > distanceP2 )
232 : 0 : p2 = i;
233 : : else
234 : 0 : p1 = i;
235 : 0 : }
236 : :
237 : : // also record count of duplicate points
238 : 0 : if ( qgsDoubleNear( distanceP1, 0 ) || qgsDoubleNear( distanceP2, 0 ) )
239 : 0 : duplicateCount++;
240 : 0 : }
241 : :
242 : 0 : if ( feedback && duplicateCount > 1 )
243 : : {
244 : 0 : feedback->pushInfo( QObject::tr( "There are at least %1 duplicate inputs, the number of output clusters may be less than was requested" ).arg( duplicateCount ) );
245 : 0 : }
246 : :
247 : : // By now two points should be found and be not the same
248 : : // Q_ASSERT( p1 != p2 && maxDistance >= 0 );
249 : :
250 : : // Accept these two points as our initial centers
251 : 0 : centers[0] = points[p1].point;
252 : 0 : centers[1] = points[p2].point;
253 : :
254 : 0 : if ( k > 2 )
255 : : {
256 : : // array of minimum distance to a point from accepted cluster centers
257 : 0 : std::vector< double > distances( n );
258 : :
259 : : // initialize array with distance to first object
260 : 0 : for ( std::size_t j = 0; j < n; j++ )
261 : : {
262 : 0 : distances[j] = points[j].point.sqrDist( centers[0] );
263 : 0 : }
264 : 0 : distances[p1] = -1;
265 : 0 : distances[p2] = -1;
266 : :
267 : : // loop i on clusters, skip 0 and 1 as found already
268 : 0 : for ( int i = 2; i < k; i++ )
269 : : {
270 : 0 : std::size_t candidateCenter = 0;
271 : 0 : double maxDistance = std::numeric_limits<double>::lowest();
272 : :
273 : : // loop j on points
274 : 0 : for ( std::size_t j = 0; j < n; j++ )
275 : : {
276 : : // accepted clusters are already marked with distance = -1
277 : 0 : if ( distances[j] < 0 )
278 : 0 : continue;
279 : :
280 : : // update minimal distance with previously accepted cluster
281 : 0 : distances[j] = std::min( points[j].point.sqrDist( centers[i - 1] ), distances[j] );
282 : :
283 : : // greedily take a point that's farthest from any of accepted clusters
284 : 0 : if ( distances[j] > maxDistance )
285 : : {
286 : 0 : candidateCenter = j;
287 : 0 : maxDistance = distances[j];
288 : 0 : }
289 : 0 : }
290 : :
291 : : // checked earlier by counting entries on input, just in case
292 : : Q_ASSERT( maxDistance >= 0 );
293 : :
294 : : // accept candidate to centers
295 : 0 : distances[candidateCenter] = -1;
296 : : // copy the point coordinates into the initial centers array
297 : 0 : centers[i] = points[candidateCenter].point;
298 : 0 : }
299 : 0 : }
300 : 0 : }
301 : :
302 : : // ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
303 : :
304 : 0 : void QgsKMeansClusteringAlgorithm::calculateKMeans( std::vector<QgsKMeansClusteringAlgorithm::Feature> &objs, std::vector<QgsPointXY> ¢ers, int k, QgsProcessingFeedback *feedback )
305 : : {
306 : 0 : int converged = false;
307 : 0 : bool changed = false;
308 : :
309 : : // avoid reallocating weights array for every iteration
310 : 0 : std::vector< uint > weights( k );
311 : :
312 : 0 : uint i = 0;
313 : 0 : for ( i = 0; i < KMEANS_MAX_ITERATIONS && !converged; i++ )
314 : : {
315 : 0 : if ( feedback && feedback->isCanceled() )
316 : 0 : break;
317 : :
318 : 0 : findNearest( objs, centers, k, changed );
319 : 0 : updateMeans( objs, centers, weights, k );
320 : 0 : converged = !changed;
321 : 0 : }
322 : :
323 : 0 : if ( !converged && feedback )
324 : 0 : feedback->reportError( QObject::tr( "Clustering did not converge after %1 iterations" ).arg( i ) );
325 : 0 : else if ( feedback )
326 : 0 : feedback->pushInfo( QObject::tr( "Clustering converged after %1 iterations" ).arg( i ) );
327 : 0 : }
328 : :
329 : : // ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
330 : :
331 : 0 : void QgsKMeansClusteringAlgorithm::findNearest( std::vector<QgsKMeansClusteringAlgorithm::Feature> &points, const std::vector<QgsPointXY> ¢ers, const int k, bool &changed )
332 : : {
333 : 0 : changed = false;
334 : 0 : std::size_t n = points.size();
335 : 0 : for ( std::size_t i = 0; i < n; i++ )
336 : : {
337 : 0 : Feature &point = points[i];
338 : :
339 : : // Initialize with distance to first cluster
340 : 0 : double currentDistance = point.point.sqrDist( centers[0] );
341 : 0 : int currentCluster = 0;
342 : :
343 : : // Check all other cluster centers and find the nearest
344 : 0 : for ( int cluster = 1; cluster < k; cluster++ )
345 : : {
346 : 0 : const double distance = point.point.sqrDist( centers[cluster] );
347 : 0 : if ( distance < currentDistance )
348 : : {
349 : 0 : currentDistance = distance;
350 : 0 : currentCluster = cluster;
351 : 0 : }
352 : 0 : }
353 : :
354 : : // Store the nearest cluster this object is in
355 : 0 : if ( point.cluster != currentCluster )
356 : : {
357 : 0 : changed = true;
358 : 0 : point.cluster = currentCluster;
359 : 0 : }
360 : 0 : }
361 : 0 : }
362 : :
363 : : // ported from https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/lwkmeans.c
364 : :
365 : 0 : void QgsKMeansClusteringAlgorithm::updateMeans( const std::vector<Feature> &points, std::vector<QgsPointXY> ¢ers, std::vector<uint> &weights, const int k )
366 : : {
367 : 0 : uint n = points.size();
368 : 0 : std::fill( weights.begin(), weights.end(), 0 );
369 : 0 : for ( int i = 0; i < k; i++ )
370 : : {
371 : 0 : centers[i].setX( 0.0 );
372 : 0 : centers[i].setY( 0.0 );
373 : 0 : }
374 : 0 : for ( uint i = 0; i < n; i++ )
375 : : {
376 : 0 : int cluster = points[i].cluster;
377 : 0 : centers[cluster] += QgsVector( points[i].point.x(),
378 : 0 : points[i].point.y() );
379 : 0 : weights[cluster] += 1;
380 : 0 : }
381 : 0 : for ( int i = 0; i < k; i++ )
382 : : {
383 : 0 : centers[i] /= weights[i];
384 : 0 : }
385 : 0 : }
386 : :
387 : :
388 : : ///@endcond
389 : :
390 : :
|
