Python计算机视觉编程学习笔记九图像分割

图像分割

（一）图割（Graph Cut）
- 1.1 从图像创建图
- 1.2 用户交互式分割
（二）利用聚类进行分割
（三）变分法

图像分割是将一幅图像分割成有意义区域的过程。区域可以是图像的前景与背景或图像中一些单独的对象。这些区域可以利用一些诸如颜色、边界或近邻相似性等特征进行构建

（一）图割（Graph Cut）

图论中的图（graph）是由若干节点（有时也称顶点）和连接节点的边构成的集合。边可以是有向的或无向的，并且这些可能有与它们相关联的权重。

图割是将一个有向图分割成两个互不相交的集合，可以用来解决很多计算机视觉方面的问题，诸如立体深度重建、图像拼接和图像分割等计算机视觉方面的不同问题。从图像像素和像素的近邻创建一个图并引入一个能量或“代价”函数，我们有可能利用图割方法将图像分割成两个或多个区域。图割的基本思想是，相似且彼此相近的像素应该划分到同一区域。

图割 C（C 是图中所有边的集合）的“代价”函数定义为所有割的边的权重求合相加：
Ecut=∑(i,j)∈CwijEcut=\sum_{(i,j)∈C}^{}wijEcut=(i,j)∈C∑wij

wijwijwij是图中节点ii到节点jj的边（i，j)（i，j)（i，j)的权重，并且是对图割CCC所有的边进行求和。

图割图像分割的思想是用图来表示图像，并对图进行划分以使割代价EcutEcutEcut最小。在用图表示图像时，增加两个额外的节点，即源点和汇点；并仅考虑那些将源点和汇点分开的割。

寻找最小割（minimum cut 或 min cut）等同于在源点和汇点间寻找最大流（maximum flow 或 max flow）
例：

from pygraph.classes.digraph import digraph
from pygraph.algorithms.minmax import maximum_flowgr = 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 为汇点，计算最大流。打印出流和割结果：

节点的标记：0 是包含图源点的部分，1 是与汇点相连的节点。

1.1 从图像创建图

给定一个邻域结构，我们可以利用图像像素作为节点定义一个图。这里讨论最简单的像素四邻域和两个图像区域（前景和背景）情况。一个四邻域（4-neighborhood）指一个像素与其正上方、正下方、左边、右边的像素直接相连。

除了像素节点外，还需要两个特定的节点——源点和汇点，来分别代表图像的前景和背景。可以利用一个简单的模型将所有像素与源点、汇点连接起来。

步骤：

每个像素节点都有一个从源点的传入边
每个像素节点都有一个到汇点的传出边
每个像素节点都有一条传入边和传出边连接到它的近邻

为了确定边的权重，需要一个能够确定这些像素点之间，像素点与源点、汇点之间边的权重（表示那条边的最大流）的分割模型。与前面一样，像素iii和像素jjj之间的边的权重记为wijwijwij，源点到像素iii的权重记为wsiwsiwsi，像素iii到汇点的权重记为witwitwit。

在前景和背景像素（从同一图像或从其他的图像）上训练出了一个贝叶斯分类器，我们就可以为前景和背景计算概率 pF(Ii)pF(Ii)pF(Ii)和 pB(Ii)pB(Ii)pB(Ii)。这里，IiIiIi是像素iii 的颜色向量。
Wsi=PF(Ii)PB(Ii)+PF(Ii)Wsi=\frac{PF(Ii)}{PB(Ii)+PF(Ii)} Wsi=PB(Ii)+PF(Ii)PF(Ii)
Wit=PB(Ii)PF(Ii)+PB(Ii)Wit=\frac{PB(Ii)}{PF(Ii)+PB(Ii)}Wit=PF(Ii)+PB(Ii)PB(Ii)
wij=ke−∣Ii−Ij∣2/σwij=ke^{−∣Ii−Ij∣^{2/σ}} wij=ke−∣Ii−Ij∣2/σ

该模型将每个像素和前景及背景（源点和汇点）连接起来，权重等于上面归一化后的概率。wijwijwij描述了近邻间像素的相似性，相似像素权重趋近于kkk，不相似的趋近于零。参数σσσ表征了随着不相似性的增加，指数次幂衰减到零的快慢。

创建 graphcut.py 的文件，从一幅图像中创建图的函数：

from pygraph.classes.digraph import digraph
from pygraph.algorithms.minmax import maximum_flow
import bayesdef 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_labels,prob = bc.classify(vim)prob_fg = prob[0]prob_bg = prob[1]#用m*n+2个节点创建图gr = digraph()gr.add_node(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 cut_graph(gr, imsize):"""    Solve max flow of graph gr and return binarylabels of the resulting segmentation."""m, n = imsizesource = m * n  # second to last is sourcesink = m * n + 1  # last is sink# cut the graphflows, cuts = maximum_flow(gr, source, sink)# convert graph to image with labelsres = zeros(m * n)for pos, label in cuts.items()[:-2]:  # don't add source/sinkres[pos] = labelreturn res.reshape((m, n))def save_as_pdf(gr, filename, show_weights=False):from pygraph.readwrite.dot import writeimport gvdot = write(gr, weighted=show_weights)gvv = gv.readstring(dot)gv.layout(gvv, 'fdp')gv.render(gvv, 'pdf', filename)def show_labeling(im, labels):"""    Show image with foreground and background areas.labels = 1 for foreground, -1 for background, 0 otherwise."""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')xticks([])yticks([])

读取图像，依据图像的两个矩形区域估算出类概率，创建一个图：

# -*- coding: utf-8 -*-
from scipy.misc import imresize
import graphcut
import cv2
import numpy as np
from pylab import  *
from PIL import Image
im = array(Image.open('D:\\Data\\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()

实验效果：

kappa的数值为2时：

kappa的数值为5时：

分析：
变量 kappa决定了近邻像素间边的相对权重。当k越大时，其分割边缘越圆润，图像的细节也逐渐丢失了。这是由于kappa决定了邻近像素间边的相对权重，权重越高，越丢失严重。

1.2 用户交互式分割

利用一些方法可以将图割分割与用户交互结合起来。例如，用户可以在一幅图像上为前景和背景提供一些标记。另一种方法是利用边界框（bounding box）或“lasso” 工具选择一个包含前景的区域。

读入一幅图像及对应的标注信息，然后将其传递到图割分割路径中：

# -*- coding: utf-8 -*-
from scipy.misc import imresize
from PCV.tools import graphcut
from PIL import Image
from numpy import *
from pylab import *def create_msr_labels(m, lasso=False):""" Create label matrix for training fromuser annotations. """labels = zeros(im.shape[:2])# backgroundlabels[m == 0] = -1labels[m == 64] = -1# foregroundif lasso:labels[m == 255] = 1else:labels[m == 128] = 1return labels# load image and annotation map
im = array(Image.open('empire.jpg'))
m = array(Image.open('empire.bmp'))
# resize
scale = 0.1
im = imresize(im, scale, interp='bilinear')
m = imresize(m, scale, interp='nearest')
# create training labels
labels = create_msr_labels(m, False)
# build graph using annotations
g = graphcut.build_bayes_graph(im, labels, kappa=2)# cut graph
res = graphcut.cut_graph(g, im.shape[:2])
# remove parts in background
res[m == 0] = 1
res[m == 64] = 1# plot the result
figure()
imshow(res)
gray()
xticks([])
yticks([])
savefig('labelplot.pdf')

分析：
左边为经过下采样后的原始图像，中间为训练的掩膜，右边为将用RGB像素值作为特征向量进行分割后的结果。观察发现，该算法对于简单图像（一根香蕉，一个小熊玩偶）的分割效果比较好，对于人像的分割，会出现部分模块分割不出来的现象，仅从与背景颜色相似度来说明人像分割不出来，似乎有些没有说服力，因为观察发现香蕉和背景颜色相似且都分割出来了，但是发现它们的背景都是纯色，并没有干扰图像分割的树，围墙，草坪，所以考虑到背景的复杂程度对图像分割可能有一定的影响。同时运行代码发现，图像越复杂，图像分割运行时间也较长。

（二）利用聚类进行分割

基于谱图理论的归一化分割算法，它将像素相似和空间近似结合起来对图像进行分割。

该方法来自定义一个分割损失函数，该损失函数不仅考虑了组的大小而且还用划分的大小对该损失函数进行“归一化”。

Encut=Ecut∑i∈AWix+Ecut∑j∈BWjxEncut= \frac{Ecut}{\sum_{i∈A}^{}W_{ix}}+\frac{Ecut}{\sum_{j∈B}^{}W_{jx}} Encut=∑i∈AWixEcut+∑j∈BWjxEcut

A 和 B 表示两个割集，并在图中分别对A 和 B 中所有其他节点（函数进行“归一化”这里指图像像素）的权重求和相加，相加求和项称为 association。对于那些像素与其他像素具有相同连接数的图像，它是对划分大小的一种粗糙度量方式。

此外归一化分割可以通过最小化下面的优化问题而求得：
Ymin=yT(D−W)yyTDyY_{min}=\frac{y^{T}(D−W)y}{y^{T}Dy}Ymin=yTDyyT(D−W)y

向量yyy包含的是离散标记，这些离散标记满足对于bbb为常数yi∈1,−byi∈{1,−b}yi∈1,−b（即yy只可以取这两个值）的约束，VVV求和为0。可以通过松弛约束条件并让 yyy 取任意实数，该问题可以变为一个容易求解的特征分解问题。这样求解的缺点是你需要对输出设定阈值或进行聚类，使它重新成为一个离散分割,问题便成为求解拉普拉斯矩阵特征向量问题：
L=D−1/2WD−1/2L=D^{−1/2}WD^{−1/2}L=D−1/2WD−1/2

难点是如何定义像素间边的权重 wijwijwij,利用原始归一化割论文中的边的权重，而连接像素iii 和像素jjj 的边的权重：

wij=e−∣Ii−Ij∣2/σge−∣xi−xj∣2/σdwij=e^{−∣Ii−Ij∣^{2}/σg}e^{−∣xi−xj∣^{2} /σd}wij=e−∣Ii−Ij∣2/σge−∣xi−xj∣2/σd

式中第一部分用来度量像素 IiIiIi和 IjIjIj之间的像素相似性，Ii(Ij)Ii(Ij)Ii(Ij)定义为 RGB 向量或灰度值；

第二部分度量图像中 xi
xi 和 xjxj 的接近程度，xi(xj)xi(xj) 定义为每个像素的坐标矢量，缩放因子 σgσgσg和 σdσdσd 决定了相对尺度和每一部件趋近 0 的快慢。

示例：建立ncut.py 的文件：

函数用来获取图像数组，并利用输入的彩色图像 RGB 值或灰度图像的灰度值创建一个特征向量。由于边的权重包含了距离部件，对于每个像素的特征向量， meshgrid() 函数来获取x 和 y 值，然后该函数会在 N 个像素上循环，并在N×N 归一化割矩阵 W 中填充值。

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

接着实现聚类过程：采用基于特征向量图像值的 K-means 聚类算法对像素进行分组：

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,full_matrices=False)# 从前 ndim 个特征向量创建特征向量# 堆叠特征向量作为矩阵的列features = array(V[:ndim]).T# k-meansfeatures = whiten(features)centroids,distortion = kmeans(features,k)code,distance = vq(features,centroids)return code,V

处理图像：

import ncut
import cv2
import numpy as np
from pylab import  *
from PIL import Image
from scipy.misc import imresize
im = array(Image.open('D:\Python\empire.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()

处理效果：

分析：
观察分割分割效果不理想，人像分割不出来，简单图像分割也会有锯齿状，考虑到采用基于特征向量图像值的 K-means 聚类算法对像素进行分组，所以调整过多次k值，依旧不理想，说明选一个合适的k值难度也比较大，对图像进行阈值处理不会给出相同的结果，对RGB值或灰度值进行聚类也一样，其中原因是它们没有考虑像素的近邻。

（三）变分法

变分问题指优化的对象是函数时，解决这类问题的算法称为变分法，如 ROF 降噪、K-means 和 SVM 。

原理： 定义关于函数的函数F[y]F[y]F[y]，亦称泛函，给定函数yyy，输出值为FFF。熵H[x]H[x]H[x]也是泛函数的一种，它定义在概率密度函数p(x)p(x)p(x)上，可等价记为H[p]H[p]H[p]。变分法类似于函数中求极值点，即寻求某个最大化或最小化泛函数F[y]F[y]F[y]的函数y(x)y(x)y(x)。利用变分法可证明两点之间的最短路径为直线以及最大熵分布为高斯分布。

Chan-Vese分割模型:

对于待分割图像区域假定一个分片常数图像模型。这里集中关注两个区域的情形，比如说是前景和背景，这个模型也可以扩展到多区域。如果用一组曲线ΓΓ将图像分离成两个区域Ω1Ω1Ω1和Ω2Ω2Ω2，分割时通过最小化Chan-Vese模型能量函数给出的：

E(Γ)=λlength(Γ)+∫Ω1(I−c1)2dx+∫Ω2(I−c2)2dxE(Γ)=λlength(Γ)+∫Ω1(I−c1)2dx+∫Ω2(I−c2)2dxE(Γ)=λlength(Γ)+∫Ω1(I−c1)2dx+∫Ω2(I−c2)2dx

变分法最小化CV模型转变成为设置阈值的ROF降噪问题:

from numpy import *
from PIL import Image
from PCV.tools import graphcut
import ncut
from scipy.misc import imresize
from pylab import *
import rofim = array(Image.open("buling.jpg").convert('L'))
U,T = rof.denoise(im,im,tolerance=0.01)
t = 0.4  #阈值
title('原始图像')
imshow(im)
import scipy.misc
scipy.misc.imsave('result.pdf',U<t*U.max())

分析：

用户交互式分割算法不能处理好复杂背景图像的分割，可以看出，图像灰度级越复杂，ROF迭代次数越多，降噪过程越久。