lwgeom_cluster_2d_kmeans()

int* lwgeom_cluster_2d_kmeans	(	const LWGEOM **	geoms,
		int	ngeoms,
		int	k
	)

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

Definition at line 77 of file lwkmeans.c.

References kmeans_config::centers, centroid(), kmeans_config::centroid_method, kmeans_config::clusters, kmeans_config::distance_method, getPoint2d_cp(), if(), kmeans_config::k, kmeans(), KMEANS_EXCEEDED_MAX_ITERATIONS, KMEANS_OK, lwalloc(), lwerror(), lwfree(), lwgeom_as_lwpoint(), lwgeom_centroid(), lwgeom_get_type(), lwgeom_is_empty(), lwkmeans_pt_centroid(), lwkmeans_pt_closest(), lwkmeans_pt_distance(), kmeans_config::max_iterations, kmeans_config::num_objs, kmeans_config::objs, LWPOINT::point, POINTTYPE, POINT2D::x, and POINT2D::y.

Referenced by ST_ClusterKMeans(), and test_kmeans().

{
        int i;
        int num_centroids = 0;
        LWGEOM **centroids;
        POINT2D *centers_raw;
        const POINT2D *cp;
        POINT2D min = { DBL_MAX,   DBL_MAX };
        POINT2D max = { -DBL_MAX, -DBL_MAX };
        double dx, dy;
        kmeans_config config;
        kmeans_result result;
        int *seen;
        int sidx = 0;

        assert(k>0);
        assert(ngeoms>0);
        assert(geoms);

    /* Initialize our static structs */
    memset(&config, 0, sizeof(kmeans_config));
    memset(&result, 0, sizeof(kmeans_result));

        if (ngeoms<k)
        {
                lwerror("%s: number of geometries is less than the number of clusters requested", __func__);
        }

        /* We'll hold the temporary centroid objects here */
        centroids = lwalloc(sizeof(LWGEOM*) * ngeoms);
        memset(centroids, 0, sizeof(LWGEOM*) * ngeoms);

        /* The vector of cluster means. We have to allocate a */
        /* chunk of memory for these because we'll be mutating them */
        /* in the kmeans algorithm */
        centers_raw = lwalloc(sizeof(POINT2D) * k);
        memset(centers_raw, 0, sizeof(POINT2D) * k);

        /* K-means configuration setup */
        config.objs = lwalloc(sizeof(Pointer) * ngeoms);
        config.num_objs = ngeoms;
        config.clusters = lwalloc(sizeof(int) * ngeoms);
        config.centers = lwalloc(sizeof(Pointer) * k);
        config.k = k;
        config.max_iterations = 0;
        config.distance_method = lwkmeans_pt_distance;
        config.centroid_method = lwkmeans_pt_centroid;

        /* Clean the memory */
        memset(config.objs, 0, sizeof(Pointer) * ngeoms);
        memset(config.clusters, 0, sizeof(int) * ngeoms);
        memset(config.centers, 0, sizeof(Pointer) * k);

        /* Prepare the list of object pointers for K-means */
        for (i = 0; i < ngeoms; i++)
        {
                const LWGEOM *geom = geoms[i];
                LWPOINT *lwpoint;

                /* Null/empty geometries get a NULL pointer */
                if ((!geom) || lwgeom_is_empty(geom))
                {
                        config.objs[i] = NULL;
                        continue;
                }

                /* If the input is a point, use its coordinates */
                /* If its not a point, convert it to one via centroid */
                if (lwgeom_get_type(geom) != POINTTYPE)
                {
                        LWGEOM *centroid = lwgeom_centroid(geom);
                        if ((!centroid) || lwgeom_is_empty(centroid))
                        {
                                config.objs[i] = NULL;
                                continue;
                        }
                        centroids[num_centroids++] = centroid;
                        lwpoint = lwgeom_as_lwpoint(centroid);
                }
                else
                {
                        lwpoint = lwgeom_as_lwpoint(geom);
                }

                /* Store a pointer to the POINT2D we are interested in */
                cp = getPoint2d_cp(lwpoint->point, 0);
                config.objs[i] = (Pointer)cp;

                /* Since we're already here, let's calculate the extrema of the set */
                if (cp->x < min.x) min.x = cp->x;
                if (cp->y < min.y) min.y = cp->y;
                if (cp->x > max.x) max.x = cp->x;
                if (cp->y > max.y) max.y = cp->y;
        }

        /*
        * We map a uniform assignment of points in the area covered by the set
        * onto actual points in the set
        */
        dx = (max.x - min.x)/k;
        dy = (max.y - min.y)/k;
        seen = lwalloc(sizeof(int)*config.k);
        memset(seen, 0, sizeof(int)*config.k);
        for (i = 0; i < k; i++)
        {
                int closest;
                POINT2D p;
                int j;

                /* Calculate a point in the range */
                p.x = min.x + dx * (i + 0.5);
                p.y = min.y + dy * (i + 0.5);

                /* Find the data point closest to the calculated point */
                closest = lwkmeans_pt_closest(config.objs, config.num_objs, &p);

                /* If something is terrible wrong w/ data, cannot find a closest */
                if (closest < 0)
                        lwerror("unable to calculate cluster seed points, too many NULLs or empties?");

                /* Ensure we aren't already using that point as a seed */
                j = 0;
                while(j < sidx)
                {
                        if (seen[j] == closest)
                        {
                                int k, t;
                                for (k = 1; k < config.num_objs; k++)
                                {
                                        t = (closest + k) % config.num_objs;
                                        if (config.objs[t])
                                        {
                                                closest = t;
                                                break;
                                        }
                                }
                                j = 0;
                        }
                        else
                        {
                                j++;
                        }
                }
                seen[sidx++] = closest;

                /* Copy the point coordinates into the initial centers array */
                /* This is ugly, but the centers array is an array of */
                /* pointers to points, not an array of points */
                centers_raw[i] = *((POINT2D*)config.objs[closest]);
                config.centers[i] = &(centers_raw[i]);
        }

        result = kmeans(&config);

        /* Before error handling, might as well clean up all the inputs */
        lwfree(config.objs);
        lwfree(config.centers);
        lwfree(centers_raw);
        lwfree(centroids);
        lwfree(seen);

        /* Good result */
        if (result == KMEANS_OK)
                return config.clusters;

        /* Bad result, not going to need the answer */
        lwfree(config.clusters);
        if (result == KMEANS_EXCEEDED_MAX_ITERATIONS)
        {
                lwerror("%s did not converge after %d iterations", __func__, config.max_iterations);
                return NULL;
        }

        /* Unknown error */
        return NULL;
}

Here is the call graph for this function:

Here is the caller graph for this function: