前言
很久以前我还在用pascal的时候，存在pascal吧的一个知名帖子，内容是使用graph库的256种颜色尝试显示各种.bmp。那个是超级简单的，因为颜色已经给定了，直接用最小二乘法算一算就没了。
后来我了解了一些图像的二值化算法（即使用黑白两种颜色表示出灰度图像），然而这个适用范围太窄了。
前几天我看到了一些YouTube上的画质极其垃圾的视频，就是经过多次有损压缩，颜色数量很少很饱和的那种。于是我突发奇想：能不能使用最多\(k\)种（自己选的）颜色，最完美地展现原来的图像？
失败的尝试
所有的值域连续的最优化题都可以尝试退火。所以我首先试着退了个火，然而效果很差，就算使用20种以上的颜色，也只有几种可以真正发挥作用。究其原因，大多数的颜色初始随机时就不会被任何一个像素点选为最接近颜色，这些颜色的微调是不会影响当前解的优劣的，就会产生很多“平台”。于是这种算法是行不通的。我很难受。
前置知识
前天看到了一个生命游戏模拟器，可以用难以想象的速度寻找一些有趣的figure。我查看了这个程序的源码，除了使用四分树的hashlife算法之外，最有趣的就是一个叫做Z曲线的算法，用来把二维的图像压缩成一维，同时保持locality。这个算法其实可以拓展至二维以上的维度。
于是我灵光一闪想到了这么一个东西：把图像里面的所有颜色记成三维空间中的点，然后问题就变成了：使用\(k\)个球，覆盖住所有的点，让最大的半径最小。然而这肯定是个我不会做的题，因为二维版本我就不会。可是我会乱搞！我把所有的点按照(前面提到的)Z-order排个序，然后在原空间中接近的点就会依然保持接近。
然后我就可以二分答案，贪心的加入，做最小球覆盖了。正确的算法应当是每次加入新球的时候使用倍增+二分确定一段序列（因为每次朴素地在\([1,n]\)中二分复杂度是不对的），然后跑随机增量。这样做的话复杂度是\(O(nm\log{nm}\log{MAX})\)的，其中\(MAX\)是距离平方最大值。可是两个log的复杂度太大了，并且绝大多数的风景图片应该不会卡最小球覆盖吧。于是我就直接一个个把点贪心地加进当前球里，直到加不下。这个最坏情况会变成四方，然而对于几乎所有的真实图片，复杂度是\(O(nm\log{MAX})\)的。
几何部分
口头说“最小球覆盖”是很简单的，然而写到代码里就变得很烦。一个点和两个点的球心是很好求的，三个点和四个点就很烦。
首先是四个点，这个部分是一道出了名的原题，按照那道题的做法，直接把点带入球的方程中，相邻两式相减，就会变成线性的方程组，所以可以高斯消元。
三个点的部分我想了很久。这是求空间三角形的外心，看起来简单的很，然而我既不会空间直线求交点，也不会求点到直线的垂线。所以普通的求外心方法全部挂掉。我尝试使用四个点的方法去做，然而这样会有一个自由元。为了加一条限制，考虑把平面方程\(Ax+By+Cz=D\)加进去一起消。所以问题转化成了求过三点的平面方程。
为了简化问题，平移使第一个点到达原点。这样显然有\(D=0\)。考虑平面上判断共线的方法：用叉积判断三角形面积是否\(=0\)即可。于是此处我们考虑混合积，同样的混合积为0即共面。混合积的形式很优美：外面是一个点积，\(A,B\)和\(C\)分别是叉积的三项系数时平面方程左侧的形式也是一个点积，所以方程系数就是叉积的系数。
程序实现
Z-order是我直接把wiki上的python实现翻译成了C++，我具体来说并不懂这段代码的原理。还有诡异的是C++不能sort一个二维数组（当cmp是一维数组比较时）。
由于高斯消元\(3\times 3\)的矩阵常数极大，图像为16万像素的时候就要跑5s。
效果
下图的左侧为原图，右侧为使用15种颜色近似的效果。图片来源为百度图片。
这次是30种颜色近似:
