Golf the K-means algorithm
C++, 479 474 bytes



Only ~20x as much as Matlab!

Golfed

#define V vector&lt;P&gt;
#define f float
struct P{f x,y,i=0;f d(P&amp;p){return(p.x-x)*(p.x-x)+(p.y-y)*(p.y-y);}f n(P&amp;p){return i?x/=i,y/=i,d(p):x=p.x,y=p.y,0;}f a(P&amp;p){x+=p.x,y+=p.y,i++;}};P z;int l(P a,P b){return a.d(z)&lt;b.d(z);}f m(f k,V&amp;p){f s=p.size(),i,j=0,t=1;V c(k),n=c,d;for(random_shuffle(p.begin(),p.end());j&lt;k;c[j].i=j++)c[j]=p[j];for(;t;c=n,n=V(k)){for(i=0;i&lt;s;i++)d=c,z=p[i],sort(d.begin(),d.end(),l),j=d[0].i,p[i].i=j,n[j].a(p[i]);for(j=t=0;j&lt;k;j++)t+=n[j].n(c[j]);}}


The input/output to the algorithm is a set of points (struct P) with x and y; and the output is the same set, with their i's set to indicate the index of the output cluster that the point ends in. 

That extra i is also used to identify the clusters. In the main loop, the closest centroid to each point is found by sorting a copy of the current centroids by proximity to that point.

This handles degenerate cases (empty clusters) by keeping the corresponding centroids' previous position (see definition of P::n, which also returns distance-to-previous-centroid). A few chars could be saved by assuming that these will not crop up.

Ungolfed, with main

#include &lt;cstdio&gt;
#include &lt;ctime&gt;
#include &lt;cstdlib&gt;
#include &lt;vector&gt;
#include &lt;algorithm&gt;
using namespace std;

#define V vector&lt;P&gt;
#define f float
struct P{
    f x,y,i=0;
    f d(P&amp;p){return(p.x-x)*(p.x-x)+(p.y-y)*(p.y-y);} // distance squared
    f n(P&amp;p){return i?x/=i,y/=i,d(p):x=p.x,y=p.y,0;} // normalize-or-reset
    f a(P&amp;p){x+=p.x,y+=p.y,i++;}                     // add coordinates
};
P z;int l(P a,P b){return a.d(z)&lt;b.d(z);}            // closer-to-z comparator 
f m(f k,V&amp;p){
    f s=p.size(),i,j=0,t=1;V c(k),n=c,d;
    for(random_shuffle(p.begin(),p.end());j&lt;k;c[j].i=j++)
        c[j]=p[j];                                // initial random assignment
    for(;t;c=n,n=V(k)){                           
        for(i=0;i&lt;s;i++)                          // assign to clusters
            d=c,z=p[i],sort(d.begin(),d.end(),l),
            j=d[0].i,p[i].i=j,n[j].a(p[i]);       // and add those coords
        for(j=t=0;j&lt;k;j++)t+=n[j].n(c[j]);        // normalize &amp; count changes
    }        
}

int main(int argc, char **argv) {
    srand((unsigned long)time(0));

    int k;
    V p;
    sscanf(argv[1], "%d", &amp;k);
    printf("Input:\n");
    for (int i=2,j=0; i&lt;argc; i+=2, j++) {
        P n;
        sscanf(argv[i], "%f", &amp;(n.x));
        sscanf(argv[i+1], "%f", &amp;(n.y));
        p.push_back(n);
        printf("%d : %f,%f\n", j, p[j].x, p[j].y);
    }

    m(k,p);
    printf("Clusters:\n");
    for (int q=0; q&lt;k; q++) {
        printf("%d\n", q);
        for (unsigned int i=0; i&lt;p.size(); i++) {
            if (p[i].i == q) printf("\t%f,%f (%d)\n", p[i].x, p[i].y, i);
        }
    }
    return 0;
}