PostGIS  2.5.0beta1dev-r@@SVN_REVISION@@
lwkmeans.c
Go to the documentation of this file.
1 /*-------------------------------------------------------------------------
2  *
3  * Copyright (c) 2018, Darafei Praliaskouski <me@komzpa.net>
4  * Copyright (c) 2016, Paul Ramsey <pramsey@cleverelephant.ca>
5  *
6  *------------------------------------------------------------------------*/
7 
8 #include <float.h>
9 #include <math.h>
10 
11 #include "liblwgeom_internal.h"
12 
13 /*
14  * When clustering lists with NULL or EMPTY elements, they will get this as
15  * their cluster number. (All the other clusters will be non-negative)
16  */
17 #define KMEANS_NULL_CLUSTER -1
18 
19 /*
20  * If the algorithm doesn't converge within this number of iterations,
21  * it will return with a failure error code.
22  */
23 #define KMEANS_MAX_ITERATIONS 1000
24 
25 static void
26 update_r(POINT2D** objs, int* clusters, uint32_t n, POINT2D** centers, uint32_t k)
27 {
28  POINT2D* obj;
29  unsigned int i;
30  double distance, curr_distance;
31  uint32_t cluster, curr_cluster;
32 
33  for (i = 0; i < n; i++)
34  {
35  obj = objs[i];
36 
37  /* Don't try to cluster NULL objects, just add them to the "unclusterable cluster" */
38  if (!obj)
39  {
40  clusters[i] = KMEANS_NULL_CLUSTER;
41  continue;
42  }
43 
44  /* Initialize with distance to first cluster */
45  curr_distance = distance2d_sqr_pt_pt(obj, centers[0]);
46  curr_cluster = 0;
47 
48  /* Check all other cluster centers and find the nearest */
49  for (cluster = 1; cluster < k; cluster++)
50  {
51  distance = distance2d_sqr_pt_pt(obj, centers[cluster]);
52  if (distance < curr_distance)
53  {
54  curr_distance = distance;
55  curr_cluster = cluster;
56  }
57  }
58 
59  /* Store the nearest cluster this object is in */
60  clusters[i] = (int) curr_cluster;
61  }
62 }
63 
64 static void
65 update_means(POINT2D** objs, int* clusters, uint32_t n, POINT2D** centers, uint32_t* weights, uint32_t k)
66 {
67  uint32_t i;
68  int cluster;
69 
70  memset(weights, 0, sizeof(uint32_t) * k);
71  for (i = 0; i < k; i++)
72  {
73  centers[i]->x = 0.0;
74  centers[i]->y = 0.0;
75  }
76  for (i = 0; i < n; i++)
77  {
78  cluster = clusters[i];
79  if (cluster != KMEANS_NULL_CLUSTER)
80  {
81  centers[cluster]->x += objs[i]->x;
82  centers[cluster]->y += objs[i]->y;
83  weights[cluster] += 1;
84  }
85  }
86  for (i = 0; i < k; i++)
87  {
88  centers[i]->x /= weights[i];
89  centers[i]->y /= weights[i];
90  }
91 }
92 
93 static int
94 kmeans(POINT2D** objs, int* clusters, uint32_t n, POINT2D** centers, uint32_t k)
95 {
96  uint32_t i = 0;
97  int* clusters_last;
98  int converged = LW_FALSE;
99  size_t clusters_sz = sizeof(int) * n;
100  uint32_t* weights;
101 
102  weights = lwalloc(sizeof(int) * k);
103 
104  /* previous cluster state array */
105  clusters_last = lwalloc(clusters_sz);
106 
107  for (i = 0; i < KMEANS_MAX_ITERATIONS && !converged; i++)
108  {
109  LW_ON_INTERRUPT(break);
110 
111  /* store the previous state of the clustering */
112  memcpy(clusters_last, clusters, clusters_sz);
113 
114  update_r(objs, clusters, n, centers, k);
115  update_means(objs, clusters, n, centers, weights, k);
116 
117  /* if all the cluster numbers are unchanged, we are at a stable solution */
118  converged = memcmp(clusters_last, clusters, clusters_sz) == 0;
119  }
120 
121  lwfree(clusters_last);
122  lwfree(weights);
123  if (!converged)
124  lwerror("%s did not converge after %d iterations", __func__, i);
125  return converged;
126 }
127 
128 static void
129 kmeans_init(POINT2D** objs, int* clusters, uint32_t n, POINT2D** centers, POINT2D* centers_raw, uint32_t k)
130 {
131  double* distances;
132  uint32_t p1 = 0, p2 = 0;
133  uint32_t i, j;
134  uint32_t duplicate_count = 1; /* a point is a duplicate of itself */
135  double max_dst = -1;
136  double dst_p1, dst_p2;
137 
138  /* k=0, k=1: "clustering" is just input validation */
139  assert(k > 1);
140 
141  /* k >= 2: find two distant points greedily */
142  for (i = 1; i < n; i++)
143  {
144  /* skip null */
145  if (!objs[i]) continue;
146 
147  /* reinit if first element happened to be null */
148  if (!objs[p1] && !objs[p2])
149  {
150  p1 = i;
151  p2 = i;
152  continue;
153  }
154 
155  /* if we found a larger distance, replace our choice */
156  dst_p1 = distance2d_sqr_pt_pt(objs[i], objs[p1]);
157  dst_p2 = distance2d_sqr_pt_pt(objs[i], objs[p2]);
158  if ((dst_p1 > max_dst) || (dst_p2 > max_dst))
159  {
160  max_dst = fmax(dst_p1, dst_p2);
161  if (dst_p1 > dst_p2)
162  p2 = i;
163  else
164  p1 = i;
165  }
166  if ((dst_p1 == 0) || (dst_p2 == 0)) duplicate_count++;
167  }
168  if (duplicate_count > 1)
169  lwnotice(
170  "%s: there are at least %u duplicate inputs, number of output clusters may be less than you requested",
171  __func__,
172  duplicate_count);
173 
174  /* by now two points should be found and non-same */
175  assert(p1 != p2 && objs[p1] && objs[p2] && max_dst >= 0);
176 
177  /* accept these two points */
178  centers_raw[0] = *((POINT2D *)objs[p1]);
179  centers[0] = &(centers_raw[0]);
180  centers_raw[1] = *((POINT2D *)objs[p2]);
181  centers[1] = &(centers_raw[1]);
182 
183  if (k > 2)
184  {
185  /* array of minimum distance to a point from accepted cluster centers */
186  distances = lwalloc(sizeof(double) * n);
187 
188  /* initialize array with distance to first object */
189  for (j = 0; j < n; j++)
190  {
191  if (objs[j])
192  distances[j] = distance2d_sqr_pt_pt(objs[j], centers[0]);
193  else
194  distances[j] = -1;
195  }
196  distances[p1] = -1;
197  distances[p2] = -1;
198 
199  /* loop i on clusters, skip 0 and 1 as found already */
200  for (i = 2; i < k; i++)
201  {
202  uint32_t candidate_center = 0;
203  double max_distance = -DBL_MAX;
204 
205  /* loop j on objs */
206  for (j = 0; j < n; j++)
207  {
208  /* empty objs and accepted clusters are already marked with distance = -1 */
209  if (distances[j] < 0) continue;
210 
211  /* update minimal distance with previosuly accepted cluster */
212  distances[j] = fmin(distance2d_sqr_pt_pt(objs[j], centers[i - 1]), distances[j]);
213 
214  /* greedily take a point that's farthest from any of accepted clusters */
215  if (distances[j] > max_distance)
216  {
217  candidate_center = j;
218  max_distance = distances[j];
219  }
220  }
221 
222  /* Checked earlier by counting entries on input, just in case */
223  assert(max_distance >= 0);
224 
225  /* accept candidate to centers */
226  distances[candidate_center] = -1;
227  /* Copy the point coordinates into the initial centers array
228  * Centers array is an array of pointers to points, not an array of points */
229  centers_raw[i] = *((POINT2D *)objs[candidate_center]);
230  centers[i] = &(centers_raw[i]);
231  }
232  lwfree(distances);
233  }
234 }
235 
236 int*
238 {
239  uint32_t i;
240  uint32_t num_centroids = 0;
241  uint32_t num_non_empty = 0;
242  LWGEOM** centroids;
243  POINT2D* centers_raw;
244  const POINT2D* cp;
245  int result = LW_FALSE;
246 
247  /* An array of objects to be analyzed.
248  * All NULL values will be returned in the KMEANS_NULL_CLUSTER. */
249  POINT2D** objs;
250 
251  /* An array of centers for the algorithm. */
252  POINT2D** centers;
253 
254  /* Array to fill in with cluster numbers. */
255  int* clusters;
256 
257  assert(k > 0);
258  assert(n > 0);
259  assert(geoms);
260 
261  if (n < k)
262  {
263  lwerror("%s: number of geometries is less than the number of clusters requested, not all clusters will get data", __func__);
264  k = n;
265  }
266 
267  /* We'll hold the temporary centroid objects here */
268  centroids = lwalloc(sizeof(LWGEOM*) * n);
269  memset(centroids, 0, sizeof(LWGEOM*) * n);
270 
271  /* The vector of cluster means. We have to allocate a chunk of memory for these because we'll be mutating them
272  * in the kmeans algorithm */
273  centers_raw = lwalloc(sizeof(POINT2D) * k);
274  memset(centers_raw, 0, sizeof(POINT2D) * k);
275 
276  /* K-means configuration setup */
277  objs = lwalloc(sizeof(POINT2D*) * n);
278  clusters = lwalloc(sizeof(int) * n);
279  centers = lwalloc(sizeof(POINT2D*) * k);
280 
281  /* Clean the memory */
282  memset(objs, 0, sizeof(POINT2D*) * n);
283  memset(clusters, 0, sizeof(int) * n);
284  memset(centers, 0, sizeof(POINT2D*) * k);
285 
286  /* Prepare the list of object pointers for K-means */
287  for (i = 0; i < n; i++)
288  {
289  const LWGEOM* geom = geoms[i];
290  LWPOINT* lwpoint;
291 
292  /* Null/empty geometries get a NULL pointer, set earlier with memset */
293  if ((!geom) || lwgeom_is_empty(geom)) continue;
294 
295  /* If the input is a point, use its coordinates */
296  /* If its not a point, convert it to one via centroid */
297  if (lwgeom_get_type(geom) != POINTTYPE)
298  {
300  if ((!centroid)) continue;
301  if (lwgeom_is_empty(centroid))
302  {
303  lwgeom_free(centroid);
304  continue;
305  }
306  centroids[num_centroids++] = centroid;
307  lwpoint = lwgeom_as_lwpoint(centroid);
308  }
309  else
310  lwpoint = lwgeom_as_lwpoint(geom);
311 
312  /* Store a pointer to the POINT2D we are interested in */
313  cp = getPoint2d_cp(lwpoint->point, 0);
314  objs[i] = (POINT2D*)cp;
315  num_non_empty++;
316  }
317 
318  if (num_non_empty < k)
319  {
320  lwnotice("%s: number of non-empty geometries is less than the number of clusters requested, not all clusters will get data", __func__);
321  k = num_non_empty;
322  }
323 
324  if (k > 1)
325  {
326  kmeans_init(objs, clusters, n, centers, centers_raw, k);
327  result = kmeans(objs, clusters, n, centers, k);
328  }
329  else
330  {
331  /* k=0: everythong is unclusterable
332  * k=1: mark up NULL and non-NULL */
333  for (i = 0; i < n; i++)
334  {
335  if (k == 0 || !objs[i])
336  clusters[i] = KMEANS_NULL_CLUSTER;
337  else
338  clusters[i] = 0;
339  }
340  result = LW_TRUE;
341  }
342 
343  /* Before error handling, might as well clean up all the inputs */
344  lwfree(objs);
345  lwfree(centers);
346  lwfree(centers_raw);
347  lwfree(centroids);
348 
349  /* Good result */
350  if (result) return clusters;
351 
352  /* Bad result, not going to need the answer */
353  lwfree(clusters);
354  return NULL;
355 }
#define KMEANS_NULL_CLUSTER
Definition: lwkmeans.c:17
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:923
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
Datum centroid(PG_FUNCTION_ARGS)
LWGEOM * lwgeom_centroid(const LWGEOM *geom)
#define LW_ON_INTERRUPT(x)
static int kmeans(POINT2D **objs, int *clusters, uint32_t n, POINT2D **centers, uint32_t k)
Definition: lwkmeans.c:94
static void update_means(POINT2D **objs, int *clusters, uint32_t n, POINT2D **centers, uint32_t *weights, uint32_t k)
Definition: lwkmeans.c:65
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
static void kmeans_init(POINT2D **objs, int *clusters, uint32_t n, POINT2D **centers, POINT2D *centers_raw, uint32_t k)
Definition: lwkmeans.c:129
POINTARRAY * point
Definition: liblwgeom.h:413
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
int * lwgeom_cluster_2d_kmeans(const LWGEOM **geoms, uint32_t n, uint32_t k)
Take a list of LWGEOMs and a number of clusters and return an integer array indicating which cluster ...
Definition: lwkmeans.c:237
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
#define KMEANS_MAX_ITERATIONS
Definition: lwkmeans.c:23
double y
Definition: liblwgeom.h:330
static void update_r(POINT2D **objs, int *clusters, uint32_t n, POINT2D **centers, uint32_t k)
Definition: lwkmeans.c:26
Datum distance(PG_FUNCTION_ARGS)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2322
void * lwalloc(size_t size)
Definition: lwutil.c:229
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1393
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364