智能剪刀——交互式图像分割
Demo Video：
注：图片来自[2]
作者：@Mapupcal
前言：
智能剪刀算法（Intelligent Scissors）是由Eric N. Mortensen和William A. Barrett于1995年提出来的一种图像分割的交互算法，可以用于2D图像分割，论文见[1][2]。
该算法可以用于辅助用户精确地勾勒出感兴趣的区域。
其可以在运行时快速地定位到该图像区域的边缘上，通过与用户交互，更精确地完成整个感兴趣区域的勾画、分割。(如上图所示)
算法理论基础：
①将图像转换成一个庞大的加权有向图。
②在运行时，通过特定的寻路算法，在出发点[a]开始，寻找一条最优路径，该路径通往用户指定的目标点[b]。
[a] 论文中为种子点，SeedPoint，为了和论文相对应，以下均使用SeedPoint，下同。
[b] 论文中为自由点，FreePoint.
注：智能剪刀算法的核心在于①，论文[1][2]着重讨论了如何确定加权有向图的边权重。
注：关于②，论文也提供了相应的算法实现，论文当中的算法就是Dijkstra最短路径搜索算法的一种变种。
将图像变成一个庞大的加权有向图：
① 图像中每个像素，在加权有向图中都将成为该图中的一个结点Node。
② 相邻的像素，互相之间存在一条加权有向边。
③ 确定每条加权有向边的权重。
注：在内存表示加权有向图的方式有很多，有每个节点基于顺序容器存放有向边的方式，也有基于矩阵的方式。在实现中，为了数据的局部性，本文采用了基于矩阵的方式去实现可提高性能。
对于某一个非边界像素p,其有八个邻居像素，如：
当转换成加权有向图之后，就变成如下这个样子：
注：上面有向边上面的数字为边的编号，并不是边的实际权重。
根据论文给出的计算边权重的计算，实际影响的因子有以下六个：
拉普拉斯交叉零点值 Fz(q)
像素点梯度值 Fg(q)
像素点梯度方向 Fd(p,q)
注：以上三个为静态因子，当不涉及Training时，有向边权重只由这三个公式确定[1][2]。
像素值 Fp(q)
“内部”像素值 Fi(q)
“外部”像素值 Fo(q)
注：以上三个为动态因子，Traning部分核心，这些因子的函数值会随用户切割图像的时候发生改变[2]，本文只会简单提及Training部分。
p表示上两图中的中心结点，q指与p相邻的结点。
[1]确定边权重的表达式如下：
其中Wz = 0.43 Wg=0.43 Wd=0.14
[2]确定边权重的表达式如下：
其中Wz = 0.3 Wg=0.3 Wd=0.1 Wp=0.1 Wi=0.1 Wo=0.1
拉普拉斯交叉零点
拉普拉斯算子可以作为图像锐化的算子，最简单的三阶3x3拉普拉斯算子如下：
用这个算子去处理图像的像素后，得到每个像素对应的拉普拉斯算子处理后的值。
对于每个结果值，我们使用如下函数去处理每一个结果值，得到拉普拉斯交叉零点的成本映射（论文中称CostMap）。
疑惑：边界像素点经过拉普拉斯算子处理后，处于边界的像素的结果值应该比较高的，而处于边界像素的有向边权重应该比较低才对。而不知为何，根据论文的公式，恰恰与这个思想相反。
注：可能是我没有找到正确的拉普拉斯算子去处理图像，我觉得论文的这个公式是反直觉的。
注：在实现中，我使用了符合我自己直觉的方法，效果还算是满意。
注：对于具备RGB彩色的图像，其像素的Fz应该为RGB三个通道中Fz的最大值。
像素的梯度值
“沿着梯度方向，函数变化得最快！”
高数的知识大多都已经丢给回高数老师了，但这句话一直都还记得。基于此，那么是不是可以认为，梯度值大的像素，处于边缘的几率就比梯度小的像素大得多呢？论文[1][2]确实也将梯度值作为权重因子考虑进去了。
梯度值计算方式个根据以下公式来计算的：
其中Ix为该像素在X方向的偏导数,Iy为该像素在Y方向的偏导数。
那么在二维图像上，如何对一个像素求其偏导数呢。
以下两个算子可用于求Ix和Iy：
注：假设横轴为x轴，竖轴为y轴，根据上面你也一定知道哪个是求x偏导，哪个是求y偏导数。（第一个为x，第二个为y）
到此为止，我们应该能够计算所有像素的对应的梯度值了。但是，这些梯度值现在还不能作为Fg的CostMap。因为对于梯度值大（有可能很大，因为求偏导的卷积运算结果是处于255*255*255~0区间）的像素点，我们要降低其梯度值对起有向边权重的影响，而对于梯度值小的，我们也要提高之。
于是我们有以下公式去达到目的：
其中
最后，让切割线段走的路径更加倾向于经过横竖像素点比斜像素点。
根据相邻结点与中心结点p的欧几里得距离，我们还要加上欧几里得距离的权重，权重如下：
注:根据欧几里得距离加权，简单来说，就是横竖边的权重应该要比斜边的权重低。（简单几何学推导就能得出该加权公式）
好了，至此，我们已完成梯度值Fg的CostMap计算。
注意：对于RBG彩色图片，请取RBG三个通道中Fg的最大值。
梯度方向
这里的梯度方向，其实也可以理解为向量PQ（对应点p、q，起点为p,终点为q）与梯度向量之间的量化权重。与梯度方向接近的边PQ权重应该越小。
注:以上均为个人猜测和理解，尚且不清楚论文的公式是否为这种意思。
见[2]中的量化公式（其中[1]此处的公式不尽相同，详情见[1]）：
其中
注：D(p)为单位向量。实现的时候记得将之归一化。
注：||p-q||表示向量的模，即向量的平方再开方。
注：对于彩色图片，在实现上，我简单使用了三通道的组合和来代替Ix,Iy。
FP、Fi、Fo
关于Fp、Fi、Fo权重因子，是和training有关的。
关于什么是training，简单点来讲，就是当用户选定了一部分分割边后（这条路径是通过智能剪刀算法来提取的），就可以用这些分割边的像素点作为Samples进行有关像素值因子的学习，从而继续影响每个像素的Fp、Fi、Fo三个CostMap的值（不仅仅是像素值相关的因子，连梯度值因子也可以做成动态的。），这样一来，每条边的权重将进行一次更新，而这次更新后的边权重，执行寻路算法后，将能更精确地选取用户想要的边界。因为有向边的权重会根据用户的交互而进行动态更新，所以说这些因子是动态的。
算法根据用户选定了的分割边的像素进行学习，这就是training。
注：动态更新每条边的权重，会增加运算操作，这些运算的成本将在每一次SeedPoint更新的时候体现。
当有向图的权重确定下来之后，就可以利用加权有向图的寻路算法来进行寻路了。
算法的思想很简单：只要每次设置SeedPoint（某个开始点）,就重新计算一次该SeedPoint到其它结点的最短路径。这个时候，所有路径的集合将构成一棵由SeedPoint扩散出去的最短路径树，有了这个最短路径树（每个结点结构中，都有一个指向路径前一个结点的指针，以便快速回溯这条最优路径），交互式的用户响应就呼之欲出了。
因为我们设置有向边权重的时候，区域边界的像素结点的边总是比其它边的权重要低，所以最短路径总是倾向于沿着区域的分割边界产生，这就得到了一把只会切割区域边界的“剪刀”。
注：注意，我用词是“倾向”而不是“总是”。具体原因你可以根据最短路径的知识来推断。
寻路算法
论文[2]中给出了寻路算法：
注：这种算法其实也算是Dijkstra最短路径搜索算法的一种变种。
注：在实现中，我使用了[3]中优化的Dijkstra算法来实现寻路算法。
快速Dijkstra算法的具体思想，就是使用一个“侵入式”的索引优先队列来从当论文[2]算法中的list，以快速获取最小成本的结点，并且，该优先队列可以非常快速地移除其底层堆中具备某个索引值的元素。
详情可以见[4]或者任何一本基础算法书当中有关最快的Dijkstra算法的实现。实现代码中使用的索引优先队列结构在[4]的211页中的习题中有所提及，个人实现是使用了[3]中提供的skeleton代码当中的优先队列和有向图结点数据结构。
实现效果
注：红色的X表示在该点设置了SeedPoint.
注：可以看到效果还是挺好的，最后来抠图。
注:关于抠图算法，在我的代码实现上，并没有对凹多边形做凸多边形分解，只是是简（偷）单（懒）地使用两点之间“扫描线”算法的实现。所以凹多边形的抠图会出现如下问题。
代码实现：智能剪刀——交互式图像分割 实现
注:项目代码基于Opencv3.2，主要用在图像卷积运算，CostMap的数据结构表示，以及图形界面的显示操作。
注：代码注重在功能的实现上，仍有较大的优化空间。
注：本代码并未实现有关training部分，如果你想知道详情，可以参考论文[2]自行实现。
Reference：
[1] Mortensen E N, Barrett W A. Intelligent scissors for image composition[C]//Proceedings of the 22nd annual conference on Computer graphics and interactive techniques. ACM, 1995: 191-198.
[2] Mortensen E N, Barrett W A. Interactive segmentation with intelligent scissors[J]. Graphical models and image processing, 1998, 60(5): 349-384.
[3] CS 4670/5670, Project 1: Image Scissors
[4] Robert Sedgewick Kevin Wayne《算法(中文第四版)》