前言
本人才疏学浅，上一篇文章中我的算法是相当复杂，效率也极低的，被打爆了。现在我们有一种迭代算法，可以在$O(k\times nmMx)$的时间复杂度内收敛到局部最优解。这里的$k$是一个（像spfa中的$k$一样的）几乎为常数的东西，$Mx$是代表点个数。它速度比我的垃圾做法快几十倍，并且效果也好很多。
参考资料
Wiki条目，请自觉翻墙。
如果您参考不了的话，我可以在此简要描述一下这个算法。
- 首先我们像之前的垃圾退火算法一样，随机确定$Mx$个点$A_1\dots A_{Mx}$，并且计算出来每一个原图中的点到最近点的距离。
- 如果$x$到$A_k$距离最近那么我们称$x$属于第$k$类。我们直接取第$i$类点的重心作为$A_i’$，即下一次迭代的$A_i$。
- 重复，直到答案变化量小于$\epsilon$，其中的$\epsilon$是一个精度有关的无穷小量。
效果
就直接拿上一篇博客中的第一个例子重现一遍吧。仍然是$15$种颜色。从上到下依次为原图、上篇博客的结果、这篇博客中的结果。
由于本人很不幸患有色弱，主观的判断也许不太靠谱。我个人认为这次的算法结果好看很多。究其原因，是最优化的目标不同。这个算法最小化的是到最近代表距离和，我上次的算法是距离最大值最小。
统计数据也是这么认为的。到最近代表距离和，图2是图3的两倍。距离最大值，图2略小于图3。
优化
有一个k-means++算法，用来求出更加优的初始值。之前的随机初始值有两个缺点：一是有某些极端情况不能收敛至最优解，二是收敛速度较慢。通过一个贪心选取两两距离尽量远的初始点，就可以大大加快收敛速度。
为了让差别更加明显，我把迭代次数从$100$改到了$10$。下图左侧为随机初始值的结果，右侧是这个贪心初始值的结果。
代码
这个代码比之前的短多了，所以直接贴上来。代码中我暴力循环了$10$次。stdin和stdout均为ppm文件。
#include<bits/stdc++.h>
using namespace std;
using ll=long long;
using uc=unsigned char;
const int k=3,Mx=10;
default_random_engine e;
int m,n,_,cnt[Mx+5],vl[1<<20],vis[1<<20];double v[Mx+5][k],sm[Mx+5][k];
struct st{int a[k];int x,y;}b[1<<20];
uc a[1010][1010][k];
template<class A,class B> double dis(A x[],B y[]){
double ans=0;
for(int i=0;i<k;i++)ans+=pow(x[i]-y[i],2);
return ans;
}
int main(){
srand(43);e.seed(rand());
scanf("P6\n%d%d\n255\n",&m,&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++){
fread(a[i][j],1,3,stdin);++_;
for(int t=0;t<3;t++)b[_].a[t]=a[i][j][t];
b[_].x=i,b[_].y=j;
}
int p=rand()%_+1;vis[p]=1;
for(int j=0;j<3;j++)v[1][j]=b[p].a[j];
for(int i=1;i<=_;i++)vl[i]=dis(v[1],b[i].a);
for(int i=2;i<=Mx;i++){
double sm=0;
for(int j=1;j<=_;j++)if(!vis[j])sm+=vl[j];
uniform_real_distribution<double> u(0,sm);
double p=u(e),now=0;
for(int j=1;j<=_;j++)if(!vis[j] && now+vl[j]>=p){
vis[j]=1;
for(int t=0;t<3;t++)v[i][t]=b[j].a[t];
for(int t=1;t<=_;t++)if(!vis[t])vl[t]=min(vl[t],(int)dis(b[t].a,v[i]));
break;
}else now+=vl[j];
}
for(int d=1;d<=10;d++){
memset(sm,0,sizeof sm);
memset(cnt,0,sizeof cnt);
for(int i=1;i<=_;i++){
double mn=1e9;vl[i]=0;
for(int j=1;j<=Mx;j++)if(dis(v[j],b[i].a)<mn)
mn=dis(v[j],b[i].a),vl[i]=j;
for(int j=0;j<k;j++)sm[vl[i]][j]+=b[i].a[j];
cnt[vl[i]]++;
}
for(int i=1;i<=Mx;i++)
for(int j=0;j<k;j++)
v[i][j]=sm[i][j]/cnt[i];
}
for(int i=1;i<=_;i++){
double mn=1e9;int mp=0;
for(int j=1;j<=Mx;j++)if(dis(v[j],b[i].a)<mn)
mn=dis(v[j],b[i].a),mp=j;
for(int j=0;j<k;j++)a[b[i].x][b[i].y][j]=v[mp][j]+.5;
}
printf("P6\n%d %d\n255\n",m,n);
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
fwrite(a[i][j],1,3,stdout);
return 0;
}
