MeanShift图像分割算法：大概是将复杂的背景，通过粗化提取整体信息，进而将图像分割。

接下来我想，将会抽出一部分时间，研究一下这个算法，以最终实现手势形状提取。

《Mean Shift: A Robust Approach Toward Feature Space Aalysis》一文中，利用Meanshift算法分割图像，大体类似于这样的效果：

看到一篇非常好哒博文

mean shift 图像分割 (一)

mean shift 图像分割（二）

mean shift 图像分割（三），讲的比较详细。

图像分割—mean shift(OpenCV源码注解)

先转发内容如下：

Reference:

[1] Mean shift: A robust approach toward feature space analysis, PAMI, 2002

[2] mean shift,非常好的ppt ，百度文库链接

[3] Pattern Recognition and Machine Learning, Bishop, 2006，Sec 2.5

[4] Computer Vision Algorithms and Applications, Richard Szeliski, 2010, Sec 5.3

[5] Kernel smoothing,MP Wand, MC Jones ,1994, Chapter 4

写在前头的话：这篇笔记看起来公式巨多，实际上只是符号表示，没啥公式推导，不过，多了就难免有差错，欢迎指正。

    Mean shitf的故事说来挺励的，早在1975年就诞生了，接着就是漫长的黑暗岁月，黑暗到几乎淡出了人们的视野，不过，命运总是善良的，95年又重新焕发生机，各种应用喷薄而出，包括目标跟踪，边缘检测，非极大值抑制等。这次就只介绍在图像分割中的应用吧，其它的我也没看。Mean shitf过程也充满正能量，描绘的是如何通过自己的努力，一步一步爬上顶峰的故事。

1 总体思想

图 1 特征空间映射：RGB图片 -> L-u特征空间

    首先meanshift是一种特征空间分析方法，要利用此方法来解决特定问题，需要将该问题映射到特征空间。对于图像分割，我们可以映射到颜色特征空间，比如将RGB图片，映射到Luv特征空间，图1是L-u二维可视化的效果。

    图像分割就是求每一个像素点的类标号。类标号取决于它在特征空间所属的cluster。对于每一个cluster，首先得有个类中心，它深深地吸引着一些点，就形成了一个类，即类中心对类中的点构成一个basin of attraction ，好比咱们的太阳系。如此，图像分割问题，就可以看成对每个像素点，找它的类中心问题，因为找到类中心就知道它是属于那一类啦，即类中心一样的点就是一类。

图2标准化后的概率密度可视化效果 -> 聚类分割结果

    密度估计的思路需要解决两个问题，what:中心是什么？how：怎么找？mean shift认为中心是概率密度（probalility density function ）的极大值点，如图2中的红色点，原文称之为mode，我这暂且用模点吧（某篇论文是如此称呼）。对于每个点怎样找到它的类中心呢？只要沿着梯度方向一步一步慢慢爬，就总能爬到极值点，图2中黑色的线，就是爬坡的轨迹。这种迭代搜索的策略在最优化中称之为 multiple restart gradient descent。不过，一般的gradient descent并不能保证收敛到局部极值，但mean shift 可以做到，因为它的步长是自适应调整的，越靠近极值点步长越小。

    也就是说meanshift的核心就两点，密度估计(Density Estimation) 和mode 搜索。对于图像数据，其分布无固定模式可循，所以密度估计必须用非参数估计，选用的是具有平滑效果的核密度估计（Kernel density estimation，KDE）。

2 算法步骤

截取这一块可视化

（a）灰度图可视化à（b）mean shift模点路径à（c）滤波后效果à（d）分割结果

    分三步走：模点搜索/图像平滑、模点聚类/合并相似区域、兼并小区域（可选）。模点搜索是为了找到每个数据点的到类中心，以中心的颜色代替自己的颜色，从而平滑图像。但模点搜索得到的模点太多，并且很多模点挨得很近，若果将每个模点都作为一类的话，类别太多，容易产生过分割，即分割太细，所以要合并掉一些模点，也就是合并相似区域。模点聚类后所得到的分割区域中，有些区域所包含的像素点太少，这些小区域也不是我们想要的，需要再次合并。

2.1 模点搜索/图像平滑

    建议先看[2]中的演示（P4-12）

    图像中的点包括两类信息：坐标空间（spatial，，），颜色空间（range ，，）。这些就构成了特征空间。

    模点搜索（OpenCV）：某一个点，它在联合特征空间中迭代搜索它的mode/模点；

    图像平滑: 将模点的颜色值赋给它自己，即.对应原文中的图像平滑，实质上是通过模点搜索，达到图像平滑的效果, 所以我合并为以一步。

    设点依次爬过的脚印为：

    出发时，它所收敛到的模点为，c代表convergence。

    第一步：如果迭代次数超过最大值（默认最多爬5次），结束搜索跳到第四步，否则，在坐标空间，筛选靠近的数据点进入下一步计算。

    OpenCV是以的坐标 为中心，边长为的方形区域内的数据点。

    其实，本应用为中心，为半径的圆形区域，那样效果更好，但是循环计算时并不方便，所以用方形区域近似。

    第二步：使用第一步幸存下来的点计算重心，并向重心移动。

    写得有点复杂了，下面解释下。是某种核函数，比如高斯分布， 是颜色空间的核平滑尺度。OpenCV使用的是最简单的均匀分布：

二维可视化效果

    是一个以（第步位置的颜色值）为球心，半径为的球体，球体内部值为1，球体外部值为0。对于经过上一步筛选后幸存的数据点，如果其颜色值满足，也就是颜色值落也在球内，那么求重心时，就要算上，否则落在球外，算重心时，就不带上它。实际上，上一步是依据坐标空间距离筛选数据点，是依据颜色距离进一步筛选数据点，上一步的筛子是矩形，这一步是球体。

    简而言之，设满足的点依次为，那么重心计算公式可以进一步化简为：

    是不是很简单呢，初中知识吧。

    注意：上文中的两个参数，是Mean shift最核心的两个参数（还有一个可选的M），具有直观的意义，分别代表坐标空间和颜色空间的核函数带宽。

    第三步：判断是否到模点了，到了就停止。

    如果，移动后颜色或者位置变化很小，则结束搜索，跳到第四步，否则重返第一步，从继续爬。

OpenCV停止搜索的条件：

    （1）坐标距离不变

    （2）颜色变化值很小。

    满足一条就可以功成身退，否则继续努力。

    第四步：将模点的颜色赋给出发点/，即。

    注意：原文这一步，不仅将模点的颜色值赋给，顺带把坐标值也赋给了，也就是说。

2.2 合并相似区域/模点聚类

    合并上一步平滑后的图像。OpenCV采用flood fill函数实现，原理很简单，看下wiki的动画就知道了，模拟洪水浸满峡谷的效果。基本上就是区域生长，从某一点出发，如果和它附近的点（4/8邻域）的颜色值相似就合并，同时再从新合并的点出发继续合并下去，直到碰到不相似的点或者该点已经属于另一类了，此时，就退回来，直到退无可退（所有的4/8邻域搜索空间都已经搜索完毕）。

    虽然很简单，但是不同的方法还是有很多需要注意的细节问题。这里假设滤波后的图像用表示。

    滤波后的两个像素点和，是否合并，可以使用颜色相似度和空间位置相似性判定。

    OpenCV只考虑颜色相似性，而忽略模点的坐标是否相似。而原算法综合了二者的信息。如果像素点，满足或者, 则这两个像素点就合并。不过OpenCV也是有考虑坐标位置的，它是只考虑原空间的4/8邻域，而原文是考虑特征空间模点的 ，相当于说OpenCV的（原空间）。

    此外，floodfill有一个特点，它不能越过已经被分类的区域，再加上没有第三步，使得OpenCV的结果，真的是惨不忍睹。原文的合并算法，具体怎么合并的还得看源代码。不过，应该不是用flood fill。

    《Computer Vision A Modern Approach》中是使用类平均距离判定是否合并。比如，，能否合并成，取决于类平均距离：

    这样做我觉得效果会更好，因为它不是单独依据边界上的两个点来判定是否合并，它是依据两个区域内部所有的点的信息综合判断。所以，它能合并两个区域，而原算法和OpenCV只能是两个点合并成一个区域，该区域又不断地合并点，一旦一个区域已经完成生长，那么它就不会和别的区域合并了。可以反证。假设先形成,区域生长的时候把给合并了，那么必定有两个点满足相似关系，连接了二者，假设这两个点为，相似，那么生长的时候就肯定已经把点合并进来了，接着把所拥有的区域全盘接收，根本不会让区域自成一类。

    当然考虑Outlier，使用中值更好。

    假设合并之后得到m类。对于原文的算法，每个像素点的标号就是其模点所属的模点集合的类标号，比如。不过，OpenCV是所属集合的类标号。

    不过，从原文结果来看，得到的结果并不是类标号，因为类标号一般都是序号，比如1，2，……，然后显示分割结果的时候，就给每一类随机分配一种独有的颜色。但原文的分割结果貌似是这一类的总体颜色值，我猜测原算法可能是用（加权）求平均的方式得到类的颜色值，然后属于这一类的像素点就用这个颜色代替。

    注意：这一步实现的是合并相似区域，但本质上还是而是合并模点，或者说模点聚类，因为每个像素点的值，就是它所属模点的颜色值/模点的联合信息。

2.3 兼并小区域

 

OpenCV的分割结果

    上一步合并了一些模点，但是，对于一些小区域，如果它和周围的颜色差异特别大，那么它们也会自成一类，这些小家伙让需要进一步合并。不过，OpenCV的实现中，并没有包含这一步，所以分割出的结果中包含了太多芝麻大点的区域，本人很不满意，有时间再加进去，还得优化下代码，这个实现实在是太慢了。怎么兼并小的区域呢？原文没说，我也没看他的源代码，我们可以直接将包含像素点少于的区域与它最相似的区域合并，实际中，小区域往往是被大区域兼并了。

3 算法原理

3.1 密度估计

    关于密度估计，这里直接使用结论，具体原理，参见第5部分:非参数密度估计。

某一点的密度估计值：

    为核函数，一般我们会使用径向对称（radially symmetric）核函数。即：

    其中为标准化常数，使得

    称为的profile,原文介绍了两种，对应两种核，这里再补充一种。

    （1）Epanechnikov Kernel

    它的profile如下：

可视化效果

    （2）Normal Kernel

    它的profile如下：

可视化效果

    （3）Uniform Kernel

    它的profile如下：

可视化效果

    3.2密度梯度估计

    3.2.1 梯度方向

    处的密度估计：

    则密度梯度估计：

    令，即这一部分又可以看成是一个核密度估计。

    物理意义：梯度方向是各个数据点的方向向量的加权求平均，即上式可以看成

蓝色圈圈—>到黄色圈圈

    例如，我们使用的是Normal Kernel，则

    想象一下几十匹马同时拉一辆车的恢宏场面，每匹马都往自己的方向拉，不过，距离越近的马，其力量越大，初中物理告诉我们，结果是合力的方向，如上图的黄色箭头。

    注意：Epanechnikov Kernel求导后实质上就是Uniform Kernel。

    3.2.2 漫漫爬坡路

    虽然，往哪个方向移动知道了，但是移动的步长并不好确定，下面转化一下形式，可以得到自适应步长:

    看起来有点复杂，实际上只是简单的替换。其中类比, 类比

    中间项的物理意义：处的核的密度估计,是求导所得，如果用Normal Kernel，则的形式和相同。

    中间项只是一个数，而最后一项就是所谓的mean shift向量，是一个方向向量，对应的就是我们的梯度方向。

    对于某一点往梯度方向移动到，则新坐标：

    物理意义：很直观，以为权值计算重心。

    当时，我们就到达了模点，由于，所以只能是。不过想要一步登天，很难，除非你出生很好，就落在模点，大多数数据点，还是得老老实实，一步一个脚印爬上去。还是设爬过的脚印依次，，则脚印公式：

    

    3.3.3 自适应步长

    可以看出步长与成反比，还是以Normal Kernel为例，越靠近模点，步长越小，反之越大。

    原文证明了，只要是凸函数，单调递减（可以不是哦），那么就能保证它总能收敛到模点，并且是单调递增的（我没看……）。只要步履不停，我们总会遇见，多么美好的世界啊，求遇见。

    3.3 图像分割领域的具体化

    本质上，mean shift解决任何问题，都是转化成密度估计问题。但具体问题还得具体分析。对于图像它有两种信息，坐标和颜色，前者为spatial 空间后者为range空间，对于单通道图片即灰度值，对于彩色图片即或者效果更好的等。二者是截然不同的属性，决定了不能等同视之。因此，我们使用多元核密度估计（multivariate kernel）。设spatial有2维，range空间，设为维。

    一元核：

    即

    图像分割中使用的多元核：

 

滤波的结果

    物理意义：分别为坐标空间核和颜色空间核的带宽（bandwidth）/尺度，我说不清，看结果吧。

    3.4回首OpenCV实现

    第二步，重心计算公式

    我们是对以为中心为边长的区域求重心，其实本应该是：

    用的是Uniform Kernel，也就是说用的是Epanechnikov Kernel

    此时，距离筛选是由核函数实现的，因此我们是对图像中所有的数据点计算重心，而不是落在为中心，为边长的区域内的点求重心。

    OpenCV的实现中， 并不是圆形的，为了循环时程序实现的方便，就用方形近似，但是严格的球体。

    不过方形的也可以写成核函数形式：

    此外，Normal Kernel 的平滑效果固然好，但是计算量大，所以主要还是用Uniform Kernel。原文说大部分场合，Uniform Kernel和Normal Kernel就能取得很好的效果。

    4延伸

    不写了，已经写得太多了……这次就只挖个坑，日后再跳

    Camshift

能够自动调整窗口的大小，能适应目标尺度变化的情况，比如人脸跟踪时，人与摄像头的距离动态变化的情况。

    带宽选择

    图像分割的带宽一般是自己调整看效果，最优带宽也能也求出来？不过，我倒想看看自适应带宽。最优带宽值看原文吧。

    Mode prune

    对于鞍点等会产生一些虚假的模点，如上图，红色线上的点可能就跑到鞍点去了，去除办法：将模点的坐标稍作移动，再从移动后的位置继续爬，如果还能爬到原来模点的位置，那就保留，否则踢掉。恩，是你的跑不了，不是你的撒手就跑。

    与双边滤波的关联

    可以看做死板的mean shift 参见[4]的5.2.1

    与分水岭分割

    逆过程，从山峰开始找山谷，参见[4]的Sec5.2.1

    补充阅读

    图像分割加速：原文提到了一种加速方法，先随机选取一部分点作为先头部队，让它们去找模点，找的过程中就会开辟出很多到模点的道路，然后呢，让其余的点插到离它最近的路走过去就好了。此外，还有层级分割的方法，OpenCV的实现应该就是其中一种实现。

    A topological approach to hierarchical segmentation using mean shift. CVPR 2007

    目标跟踪：Kernel-Based Object Tracking, PAMI 03

5 非参数密度估计

    这一部分说明为什么处的密度估计：

    其实，我觉得看bishop的那本书[2]就可以了，行云流水，精彩绝伦，其实，这本书的大部分内容都是如此精彩。我是按自己的理解写的，有些地方有改动，也会有错误，望各位看官指正。

    如果产生数据的分布形式已知，参数也已知，那么概率密度函数PDF已知，可以直接计算每一点的概率密度，比如高斯分布。如果参数不知道，那么也可以用数据估计参数，比如最小二乘估计，最大似然估计，贝叶斯参数估计等，如果连产生数据的分布形式都不知道，怎么办求概率密度呢？这就是一个非参数问题了，方法：让数据说话。

    5.1 猜一下

    对于上图中2维的情况，要估计蓝色圆域的概率密度，我相信大多数人都能凭直觉想到一种方法，那就用蓝色圈圈内的数据点个数，除以总的数据点个数，即。如果圆圈足够小，那么蓝色圈圈内部的概率密度就可以看成近似相等，那么蓝色点的概率密度应该是,是蓝色圈圈的面积。当然，也可以推广到维空间。这种算法，虽然直观，但缺乏理论支撑，下面证明，大伙的确猜对了。

    5.2理论推导

    首先说下，为什么可以用估计。

    设是一个维的数据，密度函数为，则空间中的一个区域的概率密度，即数据点落在区域的概率:

    现在假设，依据某种未知概率分布得到了N个数据点（非参数并不是无法参数化，理论上任何分布都可以参数化，毕达哥拉斯说"万物皆数"，只是参数无限维，只能当做非参数处理），则落在中的点的个数可能是,是否落在区域中就是一个二项分布：

    二项分布的期望：

    

    二项分布的方差：

    

    当时，，从参数估计角度说，前者说明是的无偏估计，后者说明是的一致估计。总之，说明，是一个很好的估计量。

    因此，

    进一步假设，比较小，那么内的可近似相等，于是：

    为的体积

    

    注意：是有偏估计，下面再说。

    由此推出，估计，有两种方法，第一种是固定看的数目，这就是kernel估计的本质（个人认为，直方图估计，Parzen windows 也是）。另外一种方法是固定看包含个数据点所需要的体积，这就是K最近邻估计。

    5.3直方图密度估计

    将数据范围划分为若干个宽度为的小栅格（bin）（也可以不等长哦），然后统计落在每个区间内的数据点个数，那么，每个区间的密度,为整个数据范围内的数据点个数。

    这个方法有很多缺陷：

    （1）第一个bin起始位置的选择会影响到结果（与bin的个数无关）

    （2）估计出来的概率密度有好多毛刺，不是连续光滑的曲线。

    （3）适合一两维的情况。维是需要的bin个数为（假设每一维都需要划分成个bin）,而且大多数bin的值为0，造成维度灾难(Curse of dimensionality)

    此外，对的大小特别敏感，小了，过拟合，不光滑，大了，太光滑，不过这是参数估计的普遍现象，前面提到的也是如此。

    5.4 K近邻密度估计（K-nearest neighbours，KNN）

    上面已经提过，，固定,看需要多大的。

    这里我们用KNN密度估计+贝叶斯 推导下KNN分类器的原理。至于怎么分类的，很简单，如果不知道的话，哈哈，看我以前写的KNN （Related部分）。

    样本属于哪一类就看它属于哪一类的可能性最大，即：

    很简单，基本的先验概率转后验概率：

    利用上面的结论，则，，

    

    所以，比较属于哪一类时，很公正，先在训练集中找K个最近的数据点，哪一类人多势众，测试样本就属于哪一类。

3类的情况

    5.5 核密度估计（kernel density estimation, KDE）

    5.5.1 Parzen windows

    点处的密度估计值,为落在以为中心的超球体的数据点个数。这与我们最开始猜测时的思想一致，只不过将超球体，换成超立方体。下面用数学符号形式化表示一下：

    好了，我们用核函数的形式表示了，这里，，为总的样本数。这种方法本质上和直方图方法没有太大的区别，Parzen windows方法是以数据点为中心，而直方图是我们自己固定好的点为中心。因此，它也会有直方图的一些缺点。比如估计的概率密度不是连续，维度灾难。

    5.5.2 Kernel smoothing

    很自然的，如果利用的数据量越大，估计出来的值就会越好，因为，我们综合的信息越多，于是我们使用所有数据点估计。采用所有样本估计的话，自然得要用加权的方法，越靠近估计点的数据点权重越大，反之，越是远离数据点，权重越小。

    前面已经介绍过具有这样属性的两种核函数。Epanechnikov Kernel和 Normal Kernel。我们可以直接替换掉，则：

    由于这两个核函数都是径向对称（radially symmetric），所以稍作了变化。

一开始，我并不理解为什么可以这么做，因为这样就已经不是窗口内的数据点个数，而是所有数据点都参合进来了，意义已经不一样了。后面我们可以通过求它的期望来进一步说明。

此外，bishop说，这个式子既可以看成，只有一个以为中心的窗口，也可以看成个以为中心的窗口，后一种介绍，我一直理解不了，但是，原文都是而不是，所以应该是第二种解释，才会这样写，我觉得第一种解释挺好，所以我都换过来了。

比如，我们使用高斯核，就有：

    注意：在靠近左边/右边的估计值有很大偏差，这是因为数据不对称，所以主要以右边/左边的数据为主，如果是回归就不会参数这种现象了。

    下面啰嗦一下上式中的

    而，所以。至于为什么要保证，下面就会知道了。

现在看看估计值的期望

    我们先做一次变量替换，。

    假设足够光滑，各阶导数都存在，我们在对在处泰勒展开一下：

    这里只推导1维的情况，维太复杂了……

    注意：无穷小项直接被我忽略了。

    第一项当然希望等于，于是我们就希望，得到第一个条件。不过，对于模式搜索来说，都可以，只要不影响到我们比较大小就好。

    第二项，等于0最好，所以我们希望

    第三项，不能发散，所以还得满足第三个条件：,原文还提到一个条件：，这个条件怎么来的，还没想清楚，很多论文也不提这个条件。

    显然这是一个有偏估计，偏差为。

    方差:

    因此，要使得期望很小，则要很小，要使方差很小则要很大。

    书上的多维推导过程，复杂，矩阵知识严重不够用。

    其中：，为了简便，我上面都是（图像分割中，一般也是如此），用来控制核函数的形状和方向，比如我们可以将高斯核改成椭圆形状。

    这里岔开一下，扯一扯目标检测。比如我们要检测图像中的椭圆形物体，用两高斯核作差，得到一个DoG（类似于墨西哥草帽），让它和图像卷积。控制它的形状和方向就能使得特定形状和方向的目标的响应值最大（和卷积核越像的区域其滤波响应值越大），从而能得到一张该目标在任何一点出现的概率图。接下来用mean shift 作模点搜索，这应该就是mean shift用于目标检测的基本原理吧，待验证。

    记录几个公式：

    ,是方阵

    ，是缩写……

    

    5.5.3直方图估计的kernel 平滑版

    参见：《Density Estimation》 Simon J. Sheather, Statistical Science 2004

    木有仔细看……

    如果假设数据服从正态分布，那么就有最优带宽，还有好多种……

    normal reference bandwidth：

    oversmoothed bandwidth

    数据的标准差，

    ：数据的个数，但是我以前看《Fast Object Detection with Entropy-Driven Evaluation》源码，用的是, 它的并不是指实际的数据，它是去掉重复后的数据，但是它论文上还是说就是样本的数目，为什么呢？

    这个用得比较多，我截取了这篇论文的部分代码，做了个小实验，……

matlab代码

close all  
ri=round (randn(500,1)*100+50);  
nb_UniQueD=numel(unique(ri));  
minScore = min(ri(:))-1;  
maxScore = max(ri(:))+1;  
scoreStd = std(ri);  sigma = 1.44 * scoreStd * nb_UniQueD^(-1/5); % not the number of sample  sigma = 1.06 * scoreStd * numel(ri)^(-1/5); % not the number of sample  numBins  = min(256,10*nb_UniQueD/2);  Sp = linspace(minScore, maxScore, numBins+1);% need to add one  H  = histc(ri, Sp);  % normalize by number of samples  Hraw = H / sum(H);  figure, subplot(211);  bar(Hraw);title('histogram estimation')  % discretization factor  discrFactor = (maxScore - minScore) / numBins;  kerSize = round(5 * sigma / discrFactor);  if kerSize(1) < 3  kerSize(1) = 3.0;  end  kerSize = double(kerSize);  % apply parzen window, kernel size such that it gets to 2 sigma  K = fspecial('gaussian', [kerSize 1], double(sigma/discrFactor) );  H = conv( Hraw, K, 'same' );  H = H + 1e-10;  H = H ./sum(H);  
subplot(212),bar(H);title('after smooth')