一个贯穿图像处理与数据挖掘的永恒问题

作者: 左飞

著有《算法之美——隐匿在数据结构背后的原理（C++版）》

原文 http://blog.csdn.net/baimafujinji/article/details/50521348

〇、序言

创新对于学术研究或产业应用都具有不言而喻的重要作用，现在国家也提出了要建立创新型国家的发展战略。如果回到我们所探讨的图像处理或数据挖掘研究，细细品读其中的某些点滴，你是否能窥探出些许启迪？首先，创新可以分成两种，一种是原始创新，另外一种就是所谓的二次创新。如果一个东西过去完全不存在，你鬼使神差的就想出来，那就是原始创新。比如图灵当初石破天惊地构想出图灵机模型就是原始创新。到现在也没有任何迹象表明，他受到了什么事或什么人的启发。事实上，现在人们（包括我学习图灵机的时候）也非常惊讶，图灵是如何提出这种前无古人的奇思妙想的！

二次创新也有很多种形式。比如逆向创新。据说人们在发明吸尘器之前最先发明的是吹尘器。一吸一吹，看似简单的一个颠倒，结果却如此神奇。现在同学们学习模式匹配算法时，必然是言必称KMP算法。的确，就原有的思路来说，KMP算法已经是做到极致了。但如果你还想有所突破呢？那就得首先打破旧有的条条框框。所以Boyer和Moore逆其道而行之，便提出了BM算法。KMP是从前向后做比较，而BM则是从后向前做比较。BM算法不仅提供了效率，更重要的是，由他们所提出的新思路为发端，后续产生了一个庞大的算法族。像BMH，BMHS等等又接踵而至。现在实际中基于BM算法的改进算法（相比于KMP）应用得其实更为广泛！

但是今天要谈的还不是逆序创新的话题。我们要谈的是二次创新中另外一种方法，暂且称之为“推衍创新”。这个思路在现代计算机科学中可谓随处可见，如果你还没有抓住他的名门，那么说明就研究工作来说，你还没入门。举一个简单的例子作为序言的结尾。最初，“伟哥”是作为治疗心绞痛的药物而研发的。但是，后来在临床试验中发现对治疗男性勃起功能障碍更加有效。所以现在它主要被应用于这方面的疾病。所以我们所说的推衍的大概意思就是，把一个领域的成果平行地转移到另外一个领域，说不定就能发挥起效！希望本文在这方面能够给大家一些启示。

一、平均值与中位数：一对死缠烂打的概念

平均数是统计学中用来衡量总体水平的一个统计量。但是，显然它并不“完美”。举个例子，现在房间里有6个人，他们的财富不尽相同，但又相差无几，这时我们可以说他们的平均身价是100万元。这个平均数基本上是有意义的，因为在假设前提下，我们知道他们6个人的财富或多或少都在100万元上下。现在马云突然来了，然后房间里变成7个人了。同样的问题，房间里所有的人平均身价可能已经突破100亿，但是这个平均数就没有什么意义了。现在房间里没有谁的身价在100亿上下。这时就引出了中位数的概念！把一组数从小到大排列，取中间位置的那个数来作为衡量该组总体水平的一个统计量。如果取包含马云在内的7个人之财富的中位数，我们知道应该还是100万左右，那么它显然是有意义的，它至少代表了这个总体中绝大多少人的情况。

你看出中位数的意义和作用了吗？现在当数据点分布比较均匀的时候，平均值是有意义的。但是一旦数据中存在异常值时，平均数就有可能失灵，这时就要用中位数来排除掉异常值的影响。但是平均数仍然有存在的价值，（只是某些时候我们要对其进行修正）。例如体育比赛时的打分机制，通常是“去掉一个最高分，去掉一个最低分，然后去平均值”。显然在体育比赛打分中，用中位数就不合适。所以我们说平均数和中位数就是一对死缠烂打的狐朋狗友！后面的内容会讨论这对概念在图像处理和数据挖掘中的重要应用。这涉及到简单平滑、中值滤波、K-means算法、k-Median算法等。你应该注意体会前面谈到的推衍创新思维。这能很好地帮助你举一反三。

二、简单平滑与中值滤波：同时联系到LeetCode上一道Hard级别的题目

现实中图像因为受到环境的影响，很容易被噪声所污染。如下图中的左上所示，这是一幅被椒盐噪声污染的图像。噪声体现为原本过渡平滑的（自然图像）区域中一个突兀的像素值。处理它最简单的策略是用一个低通滤波器对信号进行过滤。比如可以采用简单平滑算法。说白了，就是针对某个像素点，在其领域的一个小窗口内（例如3×3），对所有像素值取平均，然后用这个平均值来代替窗口中心位置的像素值。这样就能缩小噪声和非噪声像素之间的差距。也就是让原本高的值变低一点，而让原本低的值变高一点。结果如左下图所示。易见，简单平滑有一定效果，但是并不“完美”。举个例子，现在有一杯碱性溶液（PH>7），我们不断向其中加入纯水来稀释，结果PH值会越来越小。但是无论我们放多少水，这个值也不可能小于7。就算用尽全世界的水，它的整体仍然呈现碱性！

有没有更好的办法？如果你还没有想到用中位数来替代均值，那么我觉得你的头脑应该不用再继续读下去了！既然（椒盐）噪声是一个异常值，那么显然用中位数的方法来将其排掉是最好的选择了，这就是所谓的“中值”滤波的基本思想。上图右下就是采用中值滤波算法处理的图像，显然比简单平滑效果好。

但是，问题还没完！你有没有想过中值滤波相对于简单平滑的一个不足或者劣势？是的，中值滤波的复杂度太高，计算起来那是相当的耗时。为什么我会想到这个话题。因为最近我在更新我的MagicHouse（一款小型的图像处理软件http://blog.csdn.net/baimafujinji/article/details/50500757）。原来中值滤波算法是由我的刘师弟完成的。他曾经是腾讯电脑管家研发的最初核心成员，现在已经自主创业变成刘总了~我也顺便遥祝他宏图大展：）。刘同学写的中值滤波没有采用进行任何优化措施（当然这也是为了算法演示之用，如果优化了别人比较难把握原始算法的核心思想），每次执行起来都跟感觉死机了一样。于是最近我决定重写这个算法。

有兴趣的读者不妨搜一下关于中值滤波的加速算法，结论是这方面的paper很多，但我不得不告诉你大部分其实没啥创新。很多都是在串行转并行上下功夫，我真不认为这有啥意义。因为它们的基础仍然是下面我要谈的两个算法。

首先来看Leetcode上一道评级为Hard级别的题目，如下。两个有序数组，求它们合并后的中位数。这当然没啥难的，你可能会想到合并后用一个quicksort()，然后取中间位置。结果上当然可以得到正确答案。但是一定要注意：题目要求算法时间复杂度是O(log(t))，t是合并后数组的长度。但是，quicksort()的复杂度应该是O(t·log(t))。显然这样做行不通。满足时间复杂度要求才是本题的难点！

有没有什么好方法？题目给出的提示是要用“分治法”策略。而且你应该能想到是，我们要取中位数的两个子数组本来就是有序的，这个条件必须要好好利用。所以，本题的策略应该是：

该方法的核心是将原问题转变成一个寻找第k小数的问题（假设两个原序列升序排列），这样中位数实际上是第(m+n)/2小的数。所以只要解决了第k小数的问题，原问题也得以解决。 首先假设数组A和B的元素个数都大于k/2，我们比较A[k/2-1]和B[k/2-1]两个元素，这两个元素分别表示A的第k/2小的元素和B的第k/2小的元素。这两个元素比较共有三种情况：>、<和=。如果A[k/2-1]<B[k/2-1]，这表示A[0]到A[k/2-1]的元素都在A和B合并之后的前k小的元素中。换句话说，A[k/2-1]不可能大于两数组合并之后的第k小值，所以我们可以将其抛弃。反之亦然，所以当A[k/2-1]>B[k/2-1]时，我们将抛弃B[0]到B[k/2-1]的元素。

当A[k/2-1]=B[k/2-1]时，则已经找到了第k小的数，也即这个相等的元素，将其记为m。由于在A和B中分别有k/2-1个元素小于m，所以m即是第k小的数。(这里可能有人会有疑问，如果k为奇数，则m不是中位数。当然这种情况我们后面给出的代码里已有做特殊考虑，但整个算法的大体思路并无不同) 通过上面的分析，我们即可以采用递归的方式实现寻找第k小的数。此外还需要考虑几个边界条件： 如果A或者B为空，则直接返回B[k-1]或者A[k-1]； 如果k为1，我们只需要返回A[0]和B[0]中的较小值； 如果A[k/2-1]=B[k/2-1]，返回其中一个； 最终实现的代码为：

[cpp] view plain copy

1. double findKth(vector<int>& nums1, int size1, vector<int>::iterator it1,
2. vector<int>& nums2, int size2, vector<int>::iterator it2, int k)
3. {
4. //Always assume that size1 is equal or smaller than size2
5. if (size1 > size2)
6. return findKth(nums2, size2, it2, nums1, size1, it1, k);
7. if (size1 == 0) return *(it2+k-1);
9. if (k == 1) return min(*it1, *it2);
11. //Divide k into two parts
12. int offset1 = min(k / 2, size1);
13. int offset2 = k - offset1;
15. if (*(it1+offset1-1) <= *(it2+offset2-1))
16. {
17. return findKth(nums1, size1 - offset1 , it1+offset1, nums2, offset2, it2, k - offset1);
18. }
19. else
20. {
21. return findKth(nums1, offset1, it1, nums2, size2 - offset2, it2+offset2, k - offset2);
22. }
23. }
26. double findMedianSortedArrays(vector<int>& nums1, vector<int>& nums2) {
28. int m = nums1.size();
29. int n = nums2.size();
30. int total = m + n;
32. double result = findKth(nums1, m, nums1.begin(), nums2, n, nums2.begin(), (total + 1) / 2);
33. if ((total & 1) == 0) //Odd
34. {
35. result = (result + findKth( nums1, m, nums1.begin(), nums2, n, nums2.begin(), total / 2 + 1)) / 2;
36. }
37. return result;
38. }

通过对上述LeetCode题目的讨论，其实已经可以给出我们一些优化的想法了。来看下面这个图，当我们最初计算红色像素的邻域中值时，其实已经得到了红框中像素值的一个有序排列。然后在计算绿色像素的邻域中值时，我们把滑块向后移动。这时为了避免重复计算，一定要充分利用之前的有序结果。剔除最左面三个像素后的红框中的6个像素仍然有序，这时只要把新加入的绿框中的三个元素也做排序，然后得到两个有序的序列，是不是就变成了上我们讨论的问题了？而且，这个算法面对更大的滑动窗口时，优势会更为明显。

但是，如果我们只想计算3×3邻域内的中值，其实还有另外一个选择。例如下面的邻域

0 1 2

3 4 5

6 7 8

首先对窗口内的每一列分别计算最大值，中值和最小值，这样就得到了3组数据

• 最大值组：Max0 = max[P0，P3，P6]，Max1 = max[P1，P4，P7]，Max2 = max[P2，P5，P8]
• 中值组： Med0 = med[P0，P3，P6]，Med1 = med[P1，P4，P7]， Med2 = med[P2，P5，P8]
• 最小值组：Min0 = Min[P0，P3，P6]，Min1 = Min[P1，P4，P7]，Min2 = max[P2，P5，P8]

由此可以看到，最大值组中的最大值与最小值组中的最小值一定是9个元素中的最大值和最小值，不可能为中值，剩下7个；中值组中的最大值至少大于5个像素，中值组中的最小值至少小于5个像素，不可能为中值，剩下5个；最大值组中的中值至少大于5个元素，最小值组中的中值至少小于5个元素，不可能为中值，最后剩下3个要比较的元素，即 最大值组中的最小值Maxmin，中值组中的中值Medmed，最小值组中的最大值MinMax；找出这三个值中的中值为9个元素的中值。

采用上述方法，也会大大降低计算量。可见设计一个好算法永远是那么的重要！

三、数据挖掘中的K-means和K-median

聚类是将相似对象归到同一个簇中的方法，这有点像全自动分类。簇内的对象越相似，聚类的效果越好。支持向量机、神经网络所讨论的分类问题都是有监督的学习方式，现在我们所介绍的聚类则是无监督的。其中，K均值（K-means）是最基本、最简单的聚类算法。

在K均值算法中，质心是定义聚类原型（也就是机器学习获得的结果）的核心。在介绍算法实施的具体过程中，我们将演示质心的计算方法。而且你将看到除了第一次的质心是被指定的以外，此后的质心都是经由计算均值而获得的。

首先，选择K个初始质心（这K个质心并不要求来自于样本数据集），其中K是用户指定的参数，也就是所期望的簇的个数。每个数据点都被收归到距其最近之质心的分类中，而同一个质心所收归的点集为一个簇。然后，根据本次分类的结果，更新每个簇的质心。重复上述数据点分类与质心变更步骤，直到簇内数据点不再改变，或者等价地说，直到质心不再改变。

基本的K均值算法描述如下：

根据数据点到新质心的距离，再次对数据集中的数据进行分类，如图13-2(c)所示。然后，算法根据新的分类来计算新的质心，并再次根据数据点到新质心的距离，对数据集中的数据进行分类。结果发现簇内数据点不再改变，所以算法执行结束，最终的聚类结果如图13-2(d)所示。

对于距离函数和质心类型的某些组合，算法总是收敛到一个解，即K均值到达一种状态，聚类结果和质心都不再改变。但为了避免过度迭代所导致的时间消耗，实践中，也常用一个较弱的条件替换掉“质心不再发生变化”这个条件。例如，使用“直到仅有1%的点改变簇”。

尽管K均值聚类比较简单，但它也的确相当有效。它的某些变种甚至更有效， 并且不太受初始化问题的影响。但K均值并不适合所有的数据类型。它不能处理非球形簇、不同尺寸和不同密度的簇，尽管指定足够大的簇个数时它通常可以发现纯子簇。对包含离群点的数据进行聚类时，K均值也有问题。在这种情况下，离群点检测和删除大有帮助。K均值的另一个问题是，它对初值的选择是敏感的，这说明不同初值的选择所导致的迭代次数可能相差很大。此外，K值的选择也是一个问题。显然，算法本身并不能自适应地判定数据集应该被划分成几个簇。最后，K均值仅限于具有质心（均值）概念的数据。一种相关的K中心点聚类技术没有这种限制。在K中心点聚类中，我们每次选择的不再是均值，而是中位数。这种算法实现的其他细节与K均值相差不大，我们不再赘述。

最后我们给出一个实际应用的例子。（代码采用我最喜欢用做数据挖掘的R语言来实现）

一组来自世界银行的数据统计了30个国家的两项指标，我们用如下代码读入文件并显示其中最开始的几行数据。可见，数据共分散列，其中第一列是国家的名字，该项与后面的聚类分析无关，我们更关心后面两列信息。第二列给出的该国第三产业增加值占GDP的比重，最后一列给出的是人口结构中年龄大于等于65岁的人口（也就是老龄人口）占总人口的比重。

为了方便后续处理，下面对读入的数据库进行一些必要的预处理，主要是调整列标签，以及用国名替换掉行标签（同时删除包含国名的列）。

如果你绘制这些数据的散点图，不难发现这些数据大致可以分为两组。事实上，数据中有一半的国家是OECD成员国，而另外一半则属于发展中国家（包括一些东盟国家、南亚国家和拉美国家）。所以我们可以采用下面的代码来进行K均值聚类分析。

对于聚类结果，限于篇幅我们仍然只列出了最开始的几条。但是如果用图形来显示的话，可能更易于接受。下面是示例代码。

上述代码的执行结果如图13-3所示。

现在如果我问能不能提出另外一种与k-means类似的算法，你会想到什么？如果你能从k-均值算法想到提出k-中值算法，那么你算是没白读这篇文章！触类旁通，举一反三这招你算真学会了。（我想我已经无需再详细介绍k-中值算法的细节了，基本上和k-means一样，只是把所有均值出现的地方换成中值而已）这个思想看起好像很不起眼，但是你还别说，k-median算法还真的存在，而且是k-means算法的一个重要补充和改进。你可能会说提出k-median算法的人绝对是捡了一个大便宜，这么轻轻松松地就提出了一个足以留名的算法。但其实我想说，真正想到这一点的人，功力也绝对不可小觑。冰冻三尺、非一日之寒。唯有基础扎实，内力深厚的大家才能拥有这般敏锐的创新嗅觉。