文章目录
- 图像分割
- 1 图割
- 1.1 从图像创建图
- 1.2 用户交互式分割
- 2 利用聚类进行分类
- 3 变分法
图像分割
图像分割是将一幅图像分割成有意义区域的过程。区域可以是图像的前景与背景或图像中一些单独的对象。这些区域可以利用一些诸如颜色、边界或近邻相似性等特征进行构建。
1 图割
在图论中，图是由若干个节点和连接节点的边构成的集合。
图割是将一个有向图分割成两个互不相交的集合，可以用来解决诸如立体深度重建、图像拼接和图像分割等计算机视觉方面的不同问题。
图割的基本思想是：相似且彼此相近的像素应该划分到同一个区域。
图割C（C是图中所有边的集合）的“代价”函数定义为所有割的边的权重求和相加，即：
图割图像分割的思想是用图来表示图像，并对图进行划分以使割代价最小。在用图表示图像时，增加两个额外节点，即源点和汇点，并仅考虑那些将源点和汇点分开的割。
在Python计算机视觉编程中，图割使用python-graph工具包，该工具包包含许多非常有用的图算法。在cmd中键入命令安装：
pip install python-graph-core下面这个例子就是利用python-graph工具包创建的一个简单有向图并计算图的最大流/最小割：
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)得到的结果如下图所示：
图1 打印出的流和割的结果
代码首先创建有4个节点的有向图，4个节点的索引为0,1,2,3，然后用add_edge()增添边并为每条边指定特定的权重。边的权重用来衡量边的最大流容量。以节点0位源点，3位汇点，计算最大流。其中python字典包含了流穿过每条边和每个节点的标记，0是包含图源点的部分，1是与汇点相连的节点。
1.1 从图像创建图
给定一个邻域结构，可以利用图像像素作为节点定义一个图。下面给出创建图的步骤：
- 每个像素节点都有一个从源点的传入边
- 每个像素节点都有一个到汇点的传出边
- 每个像素节点都有一条传入边和传出边连接到它的近邻
为确定边的权重，需要一个能够确定这些像素点之间，像素点和源点、汇点之间边的权重的分割模型，建立模型如下：
利用这一模型，可以将每个像素和前景背景连接起来，权重等于上面归一化后的概率。
要实现这一模型，需要用到下面这个函数完成从一幅图像创建图：
def build_bayes_graph(im,labels,sigma=1e2,kappa=1):
m,n = im.shape[:2]
vim = im.reshape((-1,3))
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]
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]) + 1e-9)
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: # left exists
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: # right exists
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: # up exists
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: # down exists
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值上可以训练出一个朴素贝叶斯分类器，然后计算每一像素的分类概率，这些计算出的分类概率便是从源点出来和到汇点去的边的权重。
为了在图像可视化覆盖的标记区域，利用contourf()函数填充图像等高线间的区域，参数alpha用于设置透明度。
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')
xticks([])
yticks([])图建立起来后便需要在最优位置对图进行分割，下面这个函数可以计算最小割并将输出结果重新格式化为一个带像素标记的二值图像：
def cut_graph(gr,imsize):
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
from PIL import Image
from numpy import *
from pylab import *
from skimage.transform import resize
import numpy as np
im = array(Image.open("empire.jpg"))
"""im = imresize(im, 0.07, interp = 'bilinear')"""
im = resize(im, output_shape=(40, 56))
size = im.shape[:2]
print ("OK!!")
labels = zeros(size)
labels[3:18, 3:18] = -1
labels[-18:-3, -18:-3] = 1
print ("OK!!")
g = graphcut.build_bayes_graph(im, labels, kappa=1)
print ("OK!!")
res = graphcut.cut_graph(g, size)
print ("OK!!")
figure()
graphcut.show_labeling(im, labels)
figure()
imshow(res)
gray()
axis('off')
show()利用贝叶斯概率模型进行图像分割，最终得到的采样图像结果如下图所示：
图2 利用贝叶斯模型进行图像分割
左图为用于模型训练的标记图像，右图为分割的结果。变量kappa决定了近邻像素间边的相对权重，改变K值可以改变分割的效果。如下图所示分别为k=1，k=5的两种取值，可以看出随着K值的增大，分割边界将变得更平滑，并且细节部分也逐渐丢失。
图3 改变像素相似性和类概率之间相对权重对分割结果的影响
1.2 用户交互式分割
利用一些方法可以将图像分割与用户交互结合起来。例如用户可以在一幅图像上为前景和背景提供一些标记。
这里给出一个完整的代码示例，它会载入一幅图像及对应的标注信息，然后将其传递到我们的图像分割路径中：
from skimage.transform import resize
import graphcut
from PIL import Image
from numpy import *
from pylab import *
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('fisherman.jpg'))
m = array(Image.open('fisherman.jpg'))
im = resize(im, output_shape=(54, 38))
m = resize(m, output_shape=(54, 38))
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('1.pdf')在这段代码中，首先定义一个辅助函数用以读取这些标注图像，格式化这些标注图像便于将其传递给背景和前景训练模型，矩形框中只包含背景标记。下一步创建图并分割。由于需要有用户输入，所以要移除那些在标记背景区域里的任何前景的结果。最后绘制出分割结果。
2 利用聚类进行分类
这一部分是基于谱图理论的归一化分割算法，它将像素相似和空间近似结合起来对图像进行分割。
这一方法来自定义一个分割损失函数，该损失函数不仅需要考虑组的大小而且还用划分的大小对该损失函数进行“归一化”。A和B表示割集，并在图中分别对A和B中所有的其他节点的权重求和相加，相加求和称为association。对于那些像素与其他像素具有相同连接数的图像，它是对划分大小的一种粗糙度量方式。
定义W为边的权重矩阵，矩阵中的元素wij为连接像素i和j边的权重。D为对W每行元素求和后构成的对角矩阵，即D=diag(di)。这样，归一化分割就可以通过最小化下面的优化问题而求得：
向量y包含的是离散标记，这些离散标记满足对于b为常数yi的约束，yTD求和为0。但是由于这些约束条件，这个问题并不容易求解，需要通过松弛约束条件并让y取任意实数，则该问题变成了求解拉普拉斯矩阵特征向量问题：
将其体现在代码中：
def ncut_graph_matrix(im,sigma_d=1e2,sigma_g=1e-2):
m,n = im.shape[:2]
N = m*n
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()
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，并对这些像素进行聚类。聚类过程如下：
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))
U,sigma,V = linalg.svd(L,full_matrices=False)
features = array(V[:ndim]).T
features = whiten(features)
centroids,distortion = kmeans(features,k)
code,distance = vq(features,centroids)
return code,V接下来对这组实验进行测试：
import ncut
from skimage.transform import resize
from PIL import Image
from numpy import *
from pylab import *
im = array(Image.open('C-uniform01.ppm'))
m,n = im.shape[:2]
wid = 50
rim = resize(im, (wid, wid))
rim = array(rim, 'f')
A = ncut.ncut_graph_matrix(rim, sigma_d = 1, sigma_g = 1e-2)
code, V = ncut.cluster(A, k =10, ndim = 3)
codeim= resize(code.reshape(wid, wid), (m, n))
figure()
imshow(codeim)
gray()
show()这里用到静态手势库中的一种手势图像，并且设置聚类数k=3。得到的结果如下图所示：
图4 利用归一化分割算法分割图像
当改变不同的聚类数，得到的结果也有所不同：
图5 聚类数k=3和k=10的对比
相比较而言，聚类分割效果比用户交互式分割效果更差，细节上会有锯齿状效果，考虑到采用基于特征向量图像值的 K-means 聚类算法对像素进行分组，所以调整过多次k值，依旧不理想，说明选一个合适的k值难度也比较大，运行时间比用户交互式分割方法快一些。
上述例子中的特征向量以数组V返回，可以通过下面的代码可视化为图像：
imshow(resize(V[i].reshape(wid, wid), (m,n)))它以原图像尺寸将特征向量i显示为一幅图像。值得注意的是，对图像进行阈值处理不会给出相同的结果，对RGB值或者灰度值进行聚类也一样，原因是它们都没有考虑像素的近邻。
3 变分法
在优化过程中，当优化对象是函数时，该问题称为变分问题，解决这类问题的方法是变分法。
Chan-Vese分割模型对于待分割区域假定一个分片常数图像模型，集中关注前景和背景。若用一组曲线τ将图像分离成两个区域Ω1和Ω2，分割是通过最小化Chan-Vese模型能量函数给出的：
该能量函数用来度量与内部平均灰度常数c1和平均外部灰度常数c2的偏差。这里这两个积分是对各自区域的积分，分离曲线τ的长度用以选择更加平滑的方案。
图6 分片常数Chan-Vese分割模型
由分片常数图像U=x1c1+x2c2，可以将上式重写为（x1和x2是两区域Ω1和Ω2的特征函数）：
最小化Chan-Vese模型现在转变为设定阈值的ROF降噪问题：
import rof
from PIL import Image
from numpy import *
from pylab import *
im = array(Image.open('flower32_t0.png'))
U, T = rof.denoise(im, im, tolerance = 0.001)
t = 0.4
import scipy.misc
scipy.misc.imsave('result.pdf', U < t * U.max())在该实例中调低ROF迭代终止时的容忍阈值以确保得到足够的迭代次数。
得到的结果如下图所示：
图7 利用ROF降噪最小化Chan-Vese模型的图像分割
经过对比发现，该算法应该是三个算法中最优秀的，在该算法中，可见灰度值越复杂，迭代时间就越长。同其他算法相比，用户交互式分割算法不能处理好复杂背景图像的分割，而该算法的分割效果则更好。对于运行时间，也是三个算法中用时最短，运行最快的。
