kmeans.c

/*
 * kmeans.c
 * Map2Globe: Image to globe converter.
 *
 * Copyright (c) 2024 Bryan Franklin. All rights reserved.
 */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "kmeans.h"

double kmeans_vect_dist(kmean_vector_t *vect1, kmean_vector_t *vect2)
{
    int i=0;
    double diff=0, sum=0;

    #if 0
    printf("%s vectors:\n", __FUNCTION__);
    kmeans_print_vect(vect1);
    kmeans_print_vect(vect2);
    #endif /* 0 */

    if( vect1->length != vect2->length )
        return -1;

    for(i=0; i<vect1->length; ++i)
    {
        diff = vect1->values[i] - vect2->values[i];
        sum += diff*diff;
    }
    #if 0
    printf("%s: sum=%g; dist=%g\n", __FUNCTION__, sum, sqrt(sum));
    #endif /* 0 */

    return sqrt(sum);
}

static int kmeans_assign(kmean_vector_t *vect, kmean_vector_list_t *cents)
{
    int which=-1, i=0;
    double dist=-1, min_dist=-1;

    #if 0
    printf("vect, cents[0], cents[1]:\n");
    kmeans_print_vect(vect);
    kmeans_print_vect(&cents->data[0]);
    kmeans_print_vect(&cents->data[1]);
    #endif /* 0 */

    /* loop through all centroids, find one closest to vect */
    min_dist = dist = kmeans_vect_dist(vect,&cents->data[0]);
    which = 0;
    for(i=1; i<cents->num; ++i)
    {
        dist = kmeans_vect_dist(vect,&cents->data[i]);
        if( dist < min_dist )
        {
            min_dist = dist;
            which = i;
        }
    }

    /* return centroid id for vect */
    #if 0
    printf("returning cluster id %i\n", which);
    #endif /* 0 */
    return which;
}

static double kmeans_adjust_centers(kmean_vector_list_t *list, kmean_vector_list_t *cents)
{
    int i=0, j=0, which=0;
    int *nums=NULL;
    kmean_vector_list_t sums;

    /* allocate sums an counters */
    nums = calloc(cents->num,sizeof(int));
    kmeans_new_list(&sums,cents->num,cents->data->length);

    /* loop through list */
    for(i=0; i<list->num; ++i)
    {
        which = list->data[i].which;
        nums[which]++;
        /* sum values from each vector by class */
        for(j=0; j<sums.data[which].length; ++j)
            sums.data[which].values[j] += list->data[i].values[j];
    }

    #if 0
    printf("%s sums\n", __FUNCTION__);
    kmeans_print_list(&sums);
    #endif /* 0 */

    /* divide each sum by number in class */
    double dist = 0.0;
    kmean_vector_t old_pos;
    old_pos.values = calloc(cents->data[0].length,sizeof(double));
    for(which=0; which<cents->num; ++which)
    {
        if( nums[which] > 0 )
        {
            old_pos.length = cents->data[which].length;
            old_pos.which = cents->data[which].which;

            memcpy(old_pos.values,cents->data[which].values,sizeof(double)*cents->data[which].length);
            for(j=0; j<sums.data[which].length; ++j)
            {
                cents->data[which].values[j] = sums.data[which].values[j] / nums[which];
            }
            dist += kmeans_vect_dist(&old_pos,&cents->data[which]);
        }
    }
    free(old_pos.values); old_pos.values=NULL;

    free(nums); nums=NULL;
    kmeans_free_list(&sums);

    return dist;
}

static double kmeans_update(kmean_vector_list_t *list, kmean_vector_list_t *cents)
{
    int i=0;
    int ret = 0;

    /* assign each data point in list to nearest centroid */
    for(i=0; i<list->num; ++i)
    {
        int new_which = kmeans_assign(&list->data[i],cents);
        if( list->data[i].which != new_which )
        {
            #if 0
            printf("moving sample from cluster %i to %i.\n",
                list->data[i].which, new_which);
            #endif /* 0 */
            list->data[i].which = new_which;
        }
    }
    printf("\tmoved %i samples to different clusters.\n",ret);

    /* find new centers */
    ret = kmeans_adjust_centers(list,cents);

    return ret;
}

/* Use Lloyd's algorithm to find kmeans centers */
/* http://en.wikipedia.org/wiki/Lloyd%27s_algorithm */
int kmeans_find(kmean_vector_list_t *data, kmean_vector_list_t *cents)
{
    int iterations = 1;
    double update_dist = -1;
    while( (update_dist=kmeans_update(data,cents))>cents->num )
    {
        printf("%s: iteration %i\n", __FUNCTION__, iterations);
        kmeans_print_list(cents);
        printf("update_dist = %g\n", update_dist);
        ++iterations;
    }
    printf("%s: total iterations %i\n", __FUNCTION__, iterations);
    kmeans_print_list(cents);

    return iterations;
}

int kmeans_new_list(kmean_vector_list_t *list, int num, int len)
{
    void *ptr=NULL;
    int i=0;

    if( list==NULL )
        return -1;

    ptr = calloc(num,sizeof(kmean_vector_t));
    if(ptr==NULL )
        return -1;
    list->data = ptr;

    list->num=num;
    for(i=0; i<num; ++i)
    {
        ptr = calloc(len,sizeof(double));
        if( ptr==NULL )
            return -1;
        list->data[i].values = (double*)ptr;
        list->data[i].length = len;
    }

    return 0;
}

int kmeans_free_list(kmean_vector_list_t *list)
{
    int i=0;

    for(i=0; i<list->num; ++i)
    {
        free(list->data[i].values); list->data[i].values=NULL;
    }
    free(list->data); list->data=NULL;

    return 0;
}

int kmeans_print_vect(kmean_vector_t *vect)
{
    int j=0; 

    printf("\t<");
    for(j=0; j<vect->length; ++j)
    {
        printf("%g",vect->values[j]);
        if( j != vect->length-1 )
            printf(", ");
    }
    printf(">\n");

    return 0;
}

int kmeans_print_list(kmean_vector_list_t *list)
{
    int i=0, j=0;

    printf("Vector List:\n");
    printf("\tnum=%i\n", list->num);
    for(i=0; i<list->num; ++i)
    {
        printf("\t<");
        for(j=0; j<list->data[i].length; ++j)
        {
            printf("%g",list->data[i].values[j]);
            if( j != list->data[i].length-1 )
                printf(", ");
        }
        printf(">\n");
    }

    return 0;
}