图像分割是将一幅图像分割成有意义区域的过程。区域可以是图像的前景与背景或
图像中一些单独的对象。这些区域可以利用一些诸如颜色、边界或近邻相似性等特
征进行构建。
9.1 图割(Graph Cut)
Graph cuts是一种十分有用和流行的能量优化算法,在计算机视觉领域普遍应用于前背景分割(Image segmentation)、立体视觉(stereo vision)、抠图(Image matting)等。
此类方法把图像分割问题与图的最小割(min cut)问题相关联。首先用一个无向图G=<V,E>表示要分割的图像,V和E分别是顶点(vertex)和边(edge)的集合。此处的Graph和普通的Graph稍有不同。普通的图由顶点和边构成,如果边的有方向的,这样的图被则称为有向图,否则为无向图,且边是有权值的,不同的边可以有不同的权值,分别代表不同的物理意义。而Graph Cuts图是在普通图的基础上多了2个顶点,这2个顶点分别用符号”S”和”T”表示,统称为终端顶点。其它所有的顶点都必须和这2个顶点相连形成边集合中的一部分。所以Graph Cuts中有两种顶点,也有两种边。
第一种顶点和边是:第一种普通顶点对应于图像中的每个像素。每两个邻域顶点(对应于图像中每两个邻域像素)的连接就是一条边。这种边也叫n-links。
第二种顶点和边是:除图像像素外,还有另外两个终端顶点,叫S(source:源点,取源头之意)和T(sink:汇点,取汇聚之意)。每个普通顶点和这2个终端顶点之间都有连接,组成第二种边。这种边也叫t-links。
上图就是一个图像对应的s-t图,每个像素对应图中的一个相应顶点,另外还有s和t两个顶点。上图有两种边,实线的边表示每两个邻域普通顶点连接的边n-links,虚线的边表示每个普通顶点与s和t连接的边t-links。在前后景分割中,s一般表示前景目标,t一般表示背景。
图中每条边都有一个非负的权值we,也可以理解为cost(代价或者费用)。一个cut(割)就是图中边集合E的一个子集C,那这个割的cost(表示为|C|)就是边子集C的所有边的权值的总和。
Graph Cuts中的Cuts是指这样一个边的集合,很显然这些边集合包括了上面2种边,该集合中所有边的断开会导致残留”S”和”T”图的分开,所以就称为“割”。如果一个割,它的边的所有权值之和最小,那么这个就称为最小割,也就是图割的结果。而福特-富克森定理表明,网路的最大流max flow与最小割min cut相等。所以由Boykov和Kolmogorov发明的max-flow/min-cut算法就可以用来获得s-t图的最小割。这个最小割把图的顶点划分为两个不相交的子集S和T,其中s ∈S,t∈ T和S∪T=V 。这两个子集就对应于图像的前景像素集和背景像素集,那就相当于完成了图像分割。
图像分割可以看成pixel labeling(像素标记)问题,目标(s-node)的label设为1,背景(t-node)的label设为0,这个过程可以通过最小化图割来最小化能量函数得到。那很明显,发生在目标和背景的边界处的cut就是我们想要的(相当于把图像中背景和目标连接的地方割开,那就相当于把其分割了)。同时,这时候能量也应该是最小的。假设整幅图像的标签label(每个像素的label)为L= {l1,l2, lp },其中li为0(背景)或者1(目标)。那假设图像的分割为L时,图像的能量可以表示为:
E(L)=aR(L)+B(L)
其中,R(L)为区域项(regional term),B(L)为边界项(boundary term),而a就是区域项和边界项之间的重要因子,决定它们对能量的影响大小。如果a为0,那么就只考虑边界因素,不考虑区域因素。E(L)表示的是权值,即损失函数,也叫能量函数,图割的目标就是优化能量函数使其值达到最小。
区域项:
其中Rp(lp)表示为像素p分配标签lp的惩罚,Rp(lp)能量项的权值可以通过比较像素p的灰度和给定的目标和前景的灰度直方图来获得,换句话说就是像素p属于标签lp的概率,我希望像素p分配为其概率最大的标签lp,这时候我们希望能量最小,所以一般取概率的负对数值,故t-link的权值如下:
Rp(1) = - ln Pr(Ip|’obj’); Rp(0) = -ln Pr(Ip|’bkg’)
由上面两个公式可以看到,当像素p的灰度值属于目标的概率Pr(Ip|’obj’)大于背景Pr(Ip|’bkg’),那么Rp(1)就小于Rp(0),也就是说当像素p更有可能属于目标时,将p归类为目标就会使能量R(L)小。那么,如果全部的像素都被正确划分为目标或者背景,那么这时候能量就是最小的。
边界项:
其中,p和q为邻域像素,边界平滑项主要体现分割L的边界属性,B<p,q>可以解析为像素p和q之间不连续的惩罚,一般来说如果p和q越相似(例如它们的灰度),那么B<p,q>越大,如果他们非常不同,那么B<p,q>就接近于0。换句话说,如果两邻域像素差别很小,那么它属于同一个目标或者同一背景的可能性就很大,如果他们的差别很大,那说明这两个像素很有可能处于目标和背景的边缘部分,则被分割开的可能性比较大,所以当两邻域像素差别越大,B<p,q>越小,即能量越小。
总结:目标是将一幅图像分为目标和背景两个不相交的部分,运用图分割技术来实现。首先,图由顶点和边来组成,边有权值。那需要构建一个图,这个图有两类顶点,两类边和两类权值。普通顶点由图像每个像素组成,然后每两个邻域像素之间存在一条边,它的权值由上面说的“边界平滑能量项”来决定。还有两个终端顶点s(目标)和t(背景),每个普通顶点和s都存在连接,也就是边,边的权值由“区域能量项”Rp(1)来决定,每个普通顶点和t连接的边的权值由“区域能量项”Rp(0)来决定。这样所有边的权值就可以确定了,也就是图就确定了。这时候,就可以通过min cut算法来找到最小的割,这个min cut就是权值和最小的边的集合,这些边的断开恰好可以使目标和背景被分割开,也就是min cut对应于能量的最小化。而min cut和图的max flow是等效的,故可以通过max flow算法来找到s-t图的min cut。
在图割例子中将采用 python-graph 工具包,该工具包包含了许多非常有用的图算法,你可以在 http://code.google.com/p/python-graph/ 下载该工具包并查看文档。随后的例子里,我们要用到 maximum_flow() 函数,该函数用 Edmonds-Karp 算法(http://en.wikipedia.org/wiki/Edmonds-Karp_algorithm)计算最大流 / 最小割。采用一个完全用 Python 写成工具包的好处是安装容易且兼容性良好;不足是速度较慢。不过,对于小尺寸图像,该工具包的性能足以满足我们的需求,对于较大的图像,需要一个更快的实现。
给出一个用 python-graph 工具包计算一幅较小的图 1 的最大流 / 最小割的简单例子:
from pygraph.classes.digraph import digraph
from pygraph.algorithms.minmax import maximum_flow
gr = digraph()
gr.add_nodes([0,1,2,3])
gr.add_edge((0,1), wt=4)
gr.add_edge((1,2), wt=3)
gr.add_edge((2,3), wt=5)
gr.add_edge((0,2), wt=3)
gr.add_edge((1,3), wt=4)
flows,cuts = maximum_flow(gr,0,3)
print 'flow is:', flows
print 'cut is:', cuts
控制台输出结果:
首先,创建有 4 个节点的有向图,4 个节点的索引分别 0…3,然后用 add_edge() 增添边并为每条边指定特定的权重。边的权重用来衡量边的最大流容量。以节点 0 为源点、3 为汇点,计算最大流。上面两个 python 字典包含了流穿过每条边和每个节点的标记:0 是包含图源点的部分,1 是与汇点相连的节点。
9.1.1 从图像创建图像
给定一个邻域结构,我们可以利用图像像素作为节点定义一个图。这里我们将集中讨论最简单的像素四邻域和两个图像区域(前景和背景)情况。一个四邻域(4-neighborhood)指一个像素与其正上方、正下方、左边、右边的像素直接相连。
Graph cut的3x3图像分割示意图:我们取两个种子点(就是人为的指定分别属于目标和背景的两个像素点),然后我们建立一个图,图中边的粗细表示对应权值的大小,然后找到权值和最小的边的组合,也就是(c)中的cut,即完成了图像分割的功能。
除了像素节点外,我们还需要两个特定的节点——“源”点和“汇”点,来分别代表图像的前景和背景。我们将利用一个简单的模型将所有像素与源点、汇点连接起来。
假定我们已经在前景和背景像素(从同一图像或从其他的图像)上训练出了一个贝叶斯分类器,我们就可以为前景和背景计算概率 pF(Ii) 和 pB(Ii)。这里,Ii 是像素 i 的颜色向量。
现在我们可以为边的权重建立如下模型:
利用该模型,可以将每个像素和前景及背景(源点和汇点)连接起来,权重等于上面归一化后的概率。wij 描述了近邻间像素的相似性,相似像素权重趋近于 κ,不相似的趋近于 0。参数 σ 表征了随着不相似性的增加,指数次幂衰减到 0 的快慢。创建一个名为 graphcut.py 的文件,将下面从一幅图像创建图的函数写入该文件中:
from pygraph.classes.digraph import digraph
from pygraph.algorithms.minmax import maximum_flow
import bayes
def build_bayes_graph(im,labels,sigma=1e2,kappa=2):
""" 从像素四邻域建立一个图,前景和背景(前景用 1 标记,背景用 -1 标记,
其他的用 0 标记)由 labels 决定,并用朴素贝叶斯分类器建模 """
m,n = im.shape[:2]
# 每行是一个像素的 RGB 向量
vim = im.reshape((-1,3))
# 前景和背景(RGB)
foreground = im[labels==1].reshape((-1,3))
background = im[labels==-1].reshape((-1,3))
train_data = [foreground,background]
# 训练朴素贝叶斯分类器
bc = bayes.BayesClassifier()
bc.train(train_data)
# 获取所有像素的概率
bc_lables,prob = bc.classify(vim)
prob_fg = prob[0]
prob_bg = prob[1]
# 用m*n+2 个节点创建图
gr = digraph()
gr.add_nodes(range(m*n+2))
source = m*n # 倒数第二个是源点
sink = m*n+1 # 最后一个节点是汇点
# 归一化
for i in range(vim.shape[0]):
vim[i] = vim[i] / linalg.norm(vim[i])
# 遍历所有的节点,并添加边
for i in range(m*n):
# 从源点添加边
gr.add_edge((source,i), wt=(prob_fg[i]/(prob_fg[i]+prob_bg[i])))
# 向汇点添加边
gr.add_edge((i,sink), wt=(prob_bg[i]/(prob_fg[i]+prob_bg[i])))
# 向相邻节点添加边
if i%n != 0: # 左边存在
edge_wt = kappa*exp(-1.0*sum((vim[i]-vim[i-1])**2)/sigma)
gr.add_edge((i,i-1), wt=edge_wt)
if (i+1)%n != 0: # 如果右边存在
edge_wt = kappa*exp(-1.0*sum((vim[i]-vim[i+1])**2)/sigma)
gr.add_edge((i,i+1), wt=edge_wt)
if i//n != 0: 如果上方存在
edge_wt = kappa*exp(-1.0*sum((vim[i]-vim[i-n])**2)/sigma)
gr.add_edge((i,i-n), wt=edge_wt)
if i//n != m-1: # 如果下方存在
edge_wt = kappa*exp(-1.0*sum((vim[i]-vim[i+n])**2)/sigma)
gr.add_edge((i,i+n), wt=edge_wt)
return gr
这里我们用到了用 1 标记前景训练数据、用 -1 标记背景训练数据的一幅标记图像。基于这种标记,在 RGB 值上可以训练出一个朴素贝叶斯分类器,然后计算每一像素的分类概率,这些计算出的分类概率便是从源点出来和到汇点去的边的权重;由此可以创建一个节点为 n×m+2 的图。注意源点和汇点的索引;为了简化像素的索
引,我们将最后的两个索引作为源点和汇点的索引。
为了在图像上可视化覆盖的标记区域,我们可以利用 contourf() 函数填充图像(这里指带标记图像)等高线间的区域,参数 alpha 用于设置透明度。将下面函数增加到 graphcut.py 中:
def show_labeling(im,labels):
imshow(im)
contour(labels,[-0.5,0.5])
contourf(labels,[-1,-0.5],colors='b',alpha=0.25)
contourf(labels,[0.5,1],colors='r',alpha=0.25)
axis('off')
图建立起来后便需要在最优位置对图进行分割。下面这个函数可以计算最小割并将输出结果重新格式化为一个带像素标记的二值图像:
def cut_graph(gr,imsize):
""" 用最大流对图 gr 进行分割,并返回分割结果的二值标记 """
m,n = imsize
source = m*n # 倒数第二个节点是源点
sink = m*n+1 # 倒数第一个是汇点
# 对图进行分割
flows,cuts = maximum_flow(gr,source,sink)
# 将图转为带有标记的图像
res = zeros(m*n)
for pos,label in cuts.items()[:-2]: # 不要添加源点 / 汇点
res[pos] = label
return res.reshape((m,n))
举例会读取一幅图像,从图像的两个矩形区域估算出类概率,然后创建一个图:
from scipy.misc import imresize
import graphcut
im = array(Image.open('empire.jpg'))
im = imresize(im,0.07,interp='bilinear')
size = im.shape[:2]
# 添加两个矩形训练区域
labels = zeros(size)
labels[3:18,3:18] = -1
labels[-18:-3,-18:-3] = 1
# 创建图
g = graphcut.build_bayes_graph(im,labels,kappa=1)
# 对图进行分割
res = graphcut.cut_graph(g,size)
figure()
graphcut.show_labeling(im,labels)
figure()
imshow(res)
gray()
axis('off')
show()
我们利用了 imresize() 函数使图像小到适合我们的 Python graph 库;在该例中将图像统一缩放到原图像尺寸的 7%。图分割后将结果和训练区域一起画出来。图中图像覆盖区域为训练区域并显示出了最终的分割结果。
变量 kappa(方程中的 K)决定了近邻像素间边的相对权重。改变 K 值分割的效果如图 所示;随着 K 值增大,分割边界将变得更平滑,并且细节部分也逐步丢失。
9.1.2 用户交互式分割
利用一些方法可以将图割分割与用户交互结合起来。例如,用户可以在一幅图像上为前景和背景提供一些标记。另一种方法是利用边界框(bounding box)或“lasso” 工具选择一个包含前景的区域。
编写代码,载入一幅图像及对应的标注信息,然后将其传递到我们的图割分割路径中:
# -*- coding: utf-8 -*-
from scipy.misc import imresize
import graphcut
def create_msr_labels(m,lasso=False):
""" 从用户的注释中创建用于训练的标记矩阵 """
labels = zeros(im.shape[:2])
# 背景
labels[m==0] = -1
labels[m==64] = -1
#前景
if lasso:
labels[m == 255] = 1
else:
labels[m == 128] = 1
return labels
# 载入图像和注释图
im = array(Image.open('D:\\Python\\3.jpg'))
m = array(Image.open('D:\\Python\\3.png'))
# 调整大小
scale = 0.1
im = imresize(im, scale, interp='bilinear')
m = imresize(m, scale, interp='nearest')
# 创建训练标记
labels = create_msr_labels(m,False)
# 用注释创建图
g = graphcut.build_bayes_graph(im, labels, kappa=1)
# 图割
res = graphcut.cut_graph(g, im.shape[:2])
# 去除背景部分
res[m==0] = 1
res[m==64] = 1
# 绘制分割结果
figure()
imshow(res)
gray()
xticks([])
yticks([])
savefig('D:\\Python\\1.pdf')
首先,我们定义一个辅助函数用以读取这些标注图像,格式化这些标注图像便于将其传递给背景和前景训练模型函数,矩形框中只包含背景标记。下一步创建图并分割。由于有用户输入,所以移除那些在标记背景区域里有任何前景的结果。最后,绘制出分割结果,并通过设置这些勾选标记到一个空列表来移去这些勾选标记。这样就可以得到一个很好的边框(否则,图像中的边界在黑白图中很难看到)。
左边为经过下采样后的原始图像,中间为训练的掩膜,右边为将用RGB像素值作为特征向量进行分割后的结果。分割效果与图像复杂度成反比,即越复杂分割效果图越差,处理所耗时越长。
9.2 利用聚类进行分割
利用聚类进行分割是基于谱图理论的归一化分割算法,它将像素相似和空间近似结合起来对图像进行分割。
该方法来自定义一个分割损失函数,该损失函数不仅考虑了组的大小而且还用划分的大小对该损失函数进行“归一化”。该归一化后的分割公式将方程(9.1)中分割损失函数修改为:
A 和 B 表示两个割集,并在图中分别对 A 和 B 中所有其他节点(函数进行“归一化”这里指图像像素)的权重求和相加,相加求和项称为 association。对于那些像素与其他像素具有相同连接数的图像,它是对划分大小的一种粗糙度量方式。该算法是针对两类分割问题衍生出来的,将在下面进行讲解。
定义 W 为边的权重矩阵,矩阵中的元素 wij 为连接像素 i 和像素 j 边的权重。D 为对W 每行元素求和后构成的对角矩阵,即 D=diag(di),d w i ij = /j 。归一化分割可以通过最小化下面的优化问题而求得:
向量 y 包含的是离散标记,y(转置)D 求和为 0。由于这些约束条件,该问题不太容易求解。
然而,通过松弛约束条件并让 y 取任意实数,该问题可以变为一个容易求解的特征分解问题。这样求解的缺点是你需要对输出设定阈值或进行聚类,使它重新成为一个离散分割。
松弛该问题后,该问题便成为求解拉普拉斯矩阵特征向量问题:
正如谱聚类情形一样,现在的难点是如何定义像素间边的权重 wij。归一化割与谱聚类有很多相似之处,并且基础理论有所重叠。
我们利用原始归一化割论文 中的边的权重,通过下式给出连接像素 i 和像素 j的边的权重:
式中第一部分度量像素 Ii 和 Ij 之间的像素相似性,Ii(Ij) 定义为 RGB 向量或灰度值;
第二部分度量图像中 xi 和 xj 的接近程度,xi(xj) 定义为每个像素的坐标矢量,缩放因子 σg 和 σd 决定了相对尺度和每一部件趋近 0 的快慢。
在代码中体现:
def ncut_graph_matrix(im,sigma_d=1e2,sigma_g=1e-2):
""" 创建用于归一化割的矩阵,其中 sigma_d 和 sigma_g 是像素距离和像素相似性的权重参数 """
m,n = im.shape[:2]
N = m*n
# 归一化,并创建 RGB 或灰度特征向量
if len(im.shape)==3:
for i in range(3):
im[:,:,i] = im[:,:,i] / im[:,:,i].max()
vim = im.reshape((-1,3))
else:
im = im / im.max()
vim = im.flatten()
# x,y 坐标用于距离计算
xx,yy = meshgrid(range(n),range(m))
x,y = xx.flatten(),yy.flatten()
# 创建边线权重矩阵
W = zeros((N,N),'f')
for i in range(N):
for j in range(i,N):
d = (x[i]-x[j])**2 + (y[i]-y[j])**2
W[i,j] = W[j,i] = exp(-1.0*sum((vim[i]-vim[j])**2)/sigma_g) * exp(-d/sigma_d)
return W
该函数获取图像数组,并利用输入的彩色图像 RGB 值或灰度图像的灰度值创建一个特征向量。由于边的权重包含了距离部件,对于每个像素的特征向量,我们利用meshgrid() 函数来获取 x 和 y 值,然后该函数会在 N 个像素上循环,并在 N×N 归一化割矩阵 W 中填充值。
我们可以顺序分割每个特征向量或获取一些特征向量对它们进行聚类来计算分割结果。这里选择第二种方法,它不需要修改任意分割数也能正常工作。将拉普拉斯矩阵进行特征分解后的前 ndim 个特征向量合并在一起构成矩阵 W,并对这些像素进行聚类。下面函数实现了该聚类过程,
from scipy.cluster.vq import *
def cluster(S,k,ndim):
""" 从相似性矩阵进行谱聚类 """
# 检查对称性
if sum(abs(S-S.T)) > 1e-10:
print 'not symmetric'
# 创建拉普拉斯矩阵
rowsum = sum(abs(S),axis=0)
D = diag(1 / sqrt(rowsum + 1e-6))
L = dot(D,dot(S,D))
# 计算 L 的特征向量
U,sigma,V = linalg.svd(L)
# 从前 ndim 个特征向量创建特征向量
# 堆叠特征向量作为矩阵的列
features = array(V[:ndim]).T
# K-means 聚类
features = whiten(features)
centroids,distortion = kmeans(features,k)
code,distance = vq(features,centroids)
return code,V
这里我们采用基于特征向量图像值的 K-means 聚类算法(细节参见 6.1 节)对像素进行分组。如果你想对该结果进行实验,可以采用任一聚类算法或分组准则。现在我们可以利用该算法在一些样本图像上进行测试。下面的脚本展示了一个完整的例子:
# -*- coding: utf-8 -*-
from PCV.tools import ncut
import numpy as np
from pylab import *
from PIL import Image
from scipy.misc import imresize
im = array(Image.open('shangdalou.jpg'))
m,n = im.shape[:2]
# 调整图像的尺寸大小为 (wid,wid)
wid = 50
rim = imresize(im,(wid,wid),interp='bilinear')
rim = array(rim,'f')
# 创建归一化割矩阵
A = ncut.ncut_graph_matrix(rim,sigma_d=1,sigma_g=1e-2)
# 聚类
code,V = ncut.cluster(A,k=3,ndim=3)
# 变换到原来的图像大小
codeim = imresize(code.reshape(wid,wid),(m,n),interp='nearest')
# 绘制分割结果
figure()
imshow(codeim)
gray()
show()
因为 Numpy 的 linanlg.svd() 函数在处理大型矩阵时的计算速度并不够快(有时对于太大的矩阵甚至会给出不准确的结果),所以这里我们重新设定图像为一固定尺寸(在该例中为 50×50),以便更快地计算特征向量。在重新设定图像大小的时候我们采用了双线性插值法;因为不想插入类标记,所以在重新调整分割结果标记图像的尺寸时我们采用近邻插值法。注意,重新调整到原图像尺寸大小后第一次利用reshape 方法将一维矩阵变换为(wid,wid)二维数组。
在上面的结果中,保持所有参数都一致,只修改了图像的内容,对图像进行阈值处理不会给出相同的结果,对RGB值或灰度值进行聚类也一样,其中原因是它们没有考虑像素的近邻。
9.3 变分法
当优化的对象是函数时,该问题称为变分问题,解决这类问题的算法称为变分法。
数学原理:
函数y(x) 对于任意给定的输入变量x ,给出输出值y ;类似地,定义关于函数的函数F[y] ,亦称泛函,给定函数y ,输出值为F 。熵H[x]也是泛函的一种,它定义在概率密度函数p(x) 上,可等价记为H[p]。
泛函中变分法类似于函数中求极值点,即寻求某个最大化或最小化泛函F[y] 的函数y(x) 。利用变分法可证明两点之间的最短路径为直线以及最大熵分布为高斯分布。
对于多元函数y(x)=y(x1,…,xD) 其泰勒展开为:
对于某个泛函F[y] ,考虑y(x) 上的微小改变ϵη(x) ,其中η(x) 为任意函数,对比上式,将x 展开到无限维,从而:
为取得极值,上式一阶条件为:
注意到η(x)选取的任意性,可得泛函梯度需要处处为零,即:
下面看一个简单而有效的变分模型。Chan-Vese 分割模型对于待分割图像区域假定一个分片常数图像模型。这里我们集中关注两个区域的情形,比如前景和背景,不过这个模型也可以拓展到多区域。
如果我们用一组曲线Γ \GammaΓ将图像分离成两个区域Ω1 和 Ω2,分割是通过最小化 Chan-Vese 模型能量函数给出的:
该能量函数用来度量与内部平均灰度常数 c1 和外部平均灰度常数 c2 的偏差。这里这两个积分是对各自区域的积分,分离曲线 Γ \GammaΓ 的长度用以选择更平滑的方案。
import rof
im = array(Image.open('D:\\Python\\t0.png').convert("L"))
U,T = rof.denoise(im,im,tolerance=0.001)
t = 0.4 # 阈值
import scipy.misc
scipy.misc.imsave('D:\\Python\\t.pdf',U < t*U.max())
为确保得到足够的迭代次数,调低ROF迭代终止时的容忍阈值。利用ROF降噪最小化CV模型,得到如上面的分割结果。可以看出,图像灰度级越复杂,ROF迭代次数越多,降噪过程越久。