lwgeom_cluster_kmeans()

int * lwgeom_cluster_kmeans	(	const LWGEOM **	geoms,
		uint32_t	n,
		uint32_t	k,
		double	max_radius
	)

Take a list of LWGEOMs and a number of clusters and return an integer array indicating which cluster each geometry is in.

Parameters

geoms	the input array of LWGEOM pointers
ngeoms	the number of elements in the array
k	the number of clusters to calculate
max_radius	maximum radius of cluster before it's split

Definition at line 321 of file lwkmeans.c.

{
        uint32_t num_non_empty = 0;
 
        assert(k > 0);
        assert(n > 0);
        assert(max_radius >= 0);
        assert(geoms);
 
        if (n < k)
        {
                lwerror(
                    "%s: number of geometries is less than the number of clusters requested, not all clusters will get data",
                    __func__);
                k = n;
        }
 
        /* An array of objects to be analyzed. */
        POINT4D *objs_dense = lwalloc(sizeof(POINT4D) * n);
 
        /* Array to mark unclusterable objects. Will be returned as KMEANS_NULL_CLUSTER. */
        uint8_t *geom_valid = lwalloc(sizeof(uint8_t) * n);
        memset(geom_valid, 0, sizeof(uint8_t) * n);
 
        /* Array to fill in with cluster numbers. */
        int *clusters = lwalloc(sizeof(int) * n);
        for (uint32_t i = 0; i < n; i++)
                clusters[i] = KMEANS_NULL_CLUSTER;
 
        /* An array of clusters centers for the algorithm. */
        POINT4D *centers = lwalloc(sizeof(POINT4D) * n);
        memset(centers, 0, sizeof(POINT4D) * n);
 
        /* An array of clusters radii for the algorithm. */
        double *radii = lwalloc(sizeof(double) * n);
        memset(radii, 0, sizeof(double) * n);
 
        /* Prepare the list of object pointers for K-means */
        for (uint32_t i = 0; i < n; i++)
        {
                const LWGEOM *geom = geoms[i];
                /* Unset M values will be 1 */
                POINT4D out = {0, 0, 0, 1};
 
                /* Null/empty geometries get geom_valid=LW_FALSE set earlier with memset */
                if ((!geom) || lwgeom_is_empty(geom))
                        continue;
 
                /* If the input is a point, use its coordinates */
                if (lwgeom_get_type(geom) == POINTTYPE)
                {
                        out.x = lwpoint_get_x(lwgeom_as_lwpoint(geom));
                        out.y = lwpoint_get_y(lwgeom_as_lwpoint(geom));
                        if (lwgeom_has_z(geom))
                                out.z = lwpoint_get_z(lwgeom_as_lwpoint(geom));
                        if (lwgeom_has_m(geom))
                        {
                                out.m = lwpoint_get_m(lwgeom_as_lwpoint(geom));
                                if (out.m <= 0)
                                        lwerror("%s has an input point geometry with weight in M less or equal to 0",
                                                __func__);
                        }
                }
                else if (!lwgeom_has_z(geom))
                {
                        /* For 2D, we can take a centroid */
                        LWGEOM *centroid = lwgeom_centroid(geom);
                        if (!centroid)
                                continue;
                        if (lwgeom_is_empty(centroid))
                        {
                                lwgeom_free(centroid);
                                continue;
                        }
                        out.x = lwpoint_get_x(lwgeom_as_lwpoint(centroid));
                        out.y = lwpoint_get_y(lwgeom_as_lwpoint(centroid));
                        lwgeom_free(centroid);
                }
                else
                {
                        /* For 3D non-point, we can have a box center */
                        const GBOX *box = lwgeom_get_bbox(geom);
                        if (!gbox_is_valid(box))
                                continue;
                        out.x = (box->xmax + box->xmin) / 2;
                        out.y = (box->ymax + box->ymin) / 2;
                        out.z = (box->zmax + box->zmin) / 2;
                }
                geom_valid[i] = LW_TRUE;
                objs_dense[num_non_empty++] = out;
        }
 
        if (num_non_empty < k)
        {
                lwnotice(
                    "%s: number of non-empty geometries (%d) is less than the number of clusters (%d) requested, not all clusters will get data",
                    __func__,
                    num_non_empty,
                    k);
                k = num_non_empty;
        }
 
        uint8_t converged = LW_TRUE;
 
        if (num_non_empty > 0)
        {
                uint32_t *clusters_dense = lwalloc(sizeof(uint32_t) * num_non_empty);
                memset(clusters_dense, 0, sizeof(uint32_t) * num_non_empty);
                uint32_t output_cluster_count = kmeans(objs_dense, clusters_dense, num_non_empty, centers, radii, k, max_radius);
 
                uint32_t d = 0;
                for (uint32_t i = 0; i < n; i++)
                        if (geom_valid[i])
                                clusters[i] = (int)clusters_dense[d++];
 
                converged = output_cluster_count > 0;
                lwfree(clusters_dense);
        }
 
        /* Before error handling, might as well clean up all the inputs */
        lwfree(objs_dense);
        lwfree(centers);
        lwfree(geom_valid);
        lwfree(radii);
 
        /* Good result */
        if (converged)
                return clusters;
 
        /* Bad result, not going to need the answer */
        lwfree(clusters);
        return NULL;
}

References centroid(), gbox_is_valid(), kmeans(), KMEANS_NULL_CLUSTER, LW_TRUE, lwalloc(), lwerror(), lwfree(), lwgeom_as_lwpoint(), lwgeom_centroid(), lwgeom_free(), lwgeom_get_bbox(), lwgeom_get_type(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_is_empty(), lwnotice(), lwpoint_get_m(), lwpoint_get_x(), lwpoint_get_y(), lwpoint_get_z(), POINT4D::m, POINTTYPE, POINT4D::x, GBOX::xmax, GBOX::xmin, POINT4D::y, GBOX::ymax, GBOX::ymin, POINT4D::z, GBOX::zmax, and GBOX::zmin.

Referenced by ST_ClusterKMeans(), and test_kmeans().

Here is the call graph for this function:

Here is the caller graph for this function: