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

〇、序言

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

二次创新也有非常多种形式。比方逆向创新。

据说人们在发明吸尘器之前最先发明的是吹尘器。一吸一吹，看似简单的一个颠倒，结果却如此奇妙。

如今同学们学习模式匹配算法时，必定是言必称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级别的题目

现实中图像由于受到环境的影响，非常easy被噪声所污染。

例如以下图中的左上所看到的。这是一幅被椒盐噪声污染的图像。

噪声体现为原本过渡平滑的（自然图像）区域中一个突兀的像素值。处理它最简单的策略是用一个低通滤波器对信号进行过滤。比方能够採用简单平滑算法。说白了，就是针对某个像素点。在其领域的一个小窗体内（比如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]，返回当中一个；

终于实现的代码为：

double findKth(vector<int>& nums1, int size1, vector<int>::iterator it1, 
			   vector<int>& nums2, int size2, vector<int>::iterator it2, int k)  
{  
    //Always assume that size1 is equal or smaller than size2  
    if (size1 > size2)  
        return findKth(nums2, size2, it2, nums1, size1, it1, k);  
    if (size1 == 0)  return *(it2+k-1); 
 
    if (k == 1)  return min(*it1, *it2);  
        
    //Divide k into two parts  
    int offset1 = min(k / 2, size1);
	int offset2 = k - offset1;  
 
    if (*(it1+offset1-1) <= *(it2+offset2-1))
    {
        return findKth(nums1, size1 - offset1 , it1+offset1, nums2, offset2, it2, k - offset1);  
    }
    else
	{
    	return findKth(nums1, offset1, it1, nums2, size2 - offset2, it2+offset2, k - offset2); 
	}
}    
    

double findMedianSortedArrays(vector<int>& nums1, vector<int>& nums2) {
        
	int m = nums1.size();      
	int n = nums2.size();  
    int total = m + n; 
        
    double result = findKth(nums1, m, nums1.begin(), nums2, n, nums2.begin(), (total + 1) / 2);
	if ((total & 1) == 0)  //Odd
	{
		result = (result + 	findKth( nums1, m, nums1.begin(), nums2, n, nums2.begin(), total / 2 + 1)) / 2;
	}
	return  result;
}

通过对上述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算法的人绝对是捡了一个大廉价，这么轻轻松松地就提出了一个足以留名的算法。但事实上我想说，真正想到这一点的人。功力也绝对不可小觑。冰冻三尺、非一日之寒。唯有基础扎实。内力深厚的大家才干拥有这般敏锐的创新嗅觉。

-------------------------------------------------------------------------------------------------------------------------------------------------------

又到了广告时间了，假设你对算法感兴趣，欢迎关注我的新书《算法之美》，或者关注本博客。

详情请參见 http://blog.csdn.net/baimafujinji/article/details/50484348