Python计算机视觉——第九章 图像分割

引言

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

9.1 图割（Graph Cut）

上图由若干节点（顶点）和连接节点的边构成的集合，边可以是有向或无向的，且可能描述相关联的权重。

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

”代价“函数：图割C（所有边的集合）中所有编的权重求和相加：
$$E_{cut}=\sum _{i,j\in C}{W}_{i,j}$$
其中，系数w是节点之间的权重。

图割的思想是用图来表示图像，用图进行划分使E最小。增加两个额外的节点，即源点和汇点，仅考虑那些将源点和汇点分开的割。

寻找最小割等同于在源、汇点之间寻找最大流。

最大流基本概念：




最大流简单算法（Residual Graph）：

流量为5，但不同路径的选取会导致最大流的不同。

Ford算法：

• 选取Residul一条路径，流过之后选取反向路径
• 接着寻找路径（可以在节点选择反向路径），直到没有通路
• 容量减去残存路径，算出最大流

算法可以准确找到最大流，但是有可能时间很慢。

在图割中选择python—graph工具包，用到maxinum_flow()函数（Edmonds-Karp算法），对于小尺度图像，性能足以满足。

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个节点的有相同，索引0到3，用add_edge()增添边并指定权重，节点0为源点，3为汇点，计算最大流，打印流和割的结果：
流量=容量-残余量
割了一条边（从1到2）

Capacity

flow：

9.1.1 从图像创建图

给定领域结构，利用像素作为节点定义一个图。我们将集中讨论最简单的像素四领域和两个图像区域（前景和背景）。

除了像素节点外，还需要两个特定的节点——“源”点和“汇”点，分别代表前景和背景，主要步骤：

• 每个像素节点都有一个从源点的传入边；

• 每个像素节点都有一个到汇点的传出边；

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

为了确定边的权重，需要一个边的权重的分割模型，源点到像素i的权重记为$w_{si}$,像素i到汇点的权重记为$w_{it}$。

假设在前景和背景像素上训练出了一个贝叶斯分类器，我们就可以为前景和背景计算概率$p_F(I_i)和p_B(I_i)$,$I_i$是像素i的颜色分量。

以边为权重建立模型：

利用此模型，可将每个像素的前景和背景（源点和汇点）连接起来，权重等于归一化后的概率，$w_{ij}$描述近邻间像素的相似性，相似权重趋近于K，不相似趋近于0。

从图像中创建图的脚本：

from pylab import *
from numpy import *

from pygraph.classes.digraph import digraph
from pygraph.algorithms.minmax import maximum_flow

import bayes

""" 
Graph Cut image segmentation using max-flow/min-cut. 
"""

def build_bayes_graph(im,labels,sigma=1e2,kappa=1):
    """    Build a graph from 4-neighborhood of pixels. 
        Foreground and background is determined from
        labels (1 for foreground, -1 for background, 0 otherwise) 
        and is modeled with naive Bayes classifiers."""
    
    m,n = im.shape[:2]
    
    # RGB vector version (one pixel per row)
    vim = im.reshape((-1,3))
    
    # RGB for foreground and background
    foreground = im[labels==1].reshape((-1,3))
    background = im[labels==-1].reshape((-1,3))    
    train_data = [foreground,background]
    
    # train naive Bayes classifier
    bc = bayes.BayesClassifier()
    bc.train(train_data)

    # get probabilities for all pixels
    bc_lables,prob = bc.classify(vim)
    prob_fg = prob[0]
    prob_bg = prob[1]
    
    # create graph with m*n+2 nodes
    gr = digraph()
    gr.add_nodes(range(m*n+2))
    
    source = m*n # second to last is source
    sink = m*n+1 # last node is sink

    # normalize
    for i in range(vim.shape[0]):
        vim[i] = vim[i] / (linalg.norm(vim[i]) + 1e-9)
    
    # go through all nodes and add edges
    for i in range(m*n):
        # add edge from source
        gr.add_edge((source,i), wt=(prob_fg[i]/(prob_fg[i]+prob_bg[i])))
        
        # add edge to sink
        gr.add_edge((i,sink), wt=(prob_bg[i]/(prob_fg[i]+prob_bg[i])))
        
        # add edges to neighbors
        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    
        
    
def cut_graph(gr,imsize):
    """    Solve max flow of graph gr and return binary 
        labels of the resulting segmentation."""
    
    m,n = imsize
    source = m*n # second to last is source
    sink = m*n+1 # last is sink
    
    # cut the graph
    flows,cuts = maximum_flow(gr,source,sink)
    
    # convert graph to image with labels
    res = zeros(m*n)
    for pos,label in list(cuts.items())[:-2]: #don't add source/sink
        res[pos] = label

    return res.reshape((m,n))
    
    
def save_as_pdf(gr,filename,show_weights=False):

    from pygraph.readwrite.dot import write
    import gv
    dot = 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([])

1标记前景训练数据，-1标记背景训练数据，在RGB值上训练出分类器，返回标记和概率，分类的概率是边的权重，由此创建节点为$n\times m+2$的图，最后两个索引是源点和汇点。

show_labeling（）用于可视化覆盖大的标记区域，利用contourf()函数填充图像等高线的区域，alpha用于设置透明度。

前景labels = 1,背景labels = -1,其他为0。

cut_graph()计算最小割并将输出结果变换为带像素标记的二值图像。cut.item[:-2]函数作用是返回字典的列表，{‘pos’:label}，忽略最后两个元素。

下面的脚本将分割图像，从矩形的两个区域估算类概率：

import cv2 as cv
import graphcut
from PIL import Image
from numpy import *


im = array(Image.open('pic/empire.jpg'))

im = cv.resize(im, None, fx = 0.07, fy = 0.07, interpolation = cv.INTER_NEAREST)


size = im.shape[:2]

# 添加两个矩形训练区域
labels = zeros(size)
# 索引至17
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()

Keppa决定了近邻像素间的相对权重，随着K值增大，边界将变得更平滑，细节部分也逐渐丢失。

9.2 利用聚类进行分割

本节，分割方法基于谱图理论的归一化分割算法，将像素相似和空间近似结合起来对图像进行分割。
该方法来自定义分割损失函数，不仅考虑组的大小，还用划分的大小对该损失函数进行归一化,归一化后：
$$E_{ncut}=\frac{E_{cut}}{{\sum }_{i\in A}{w}_{ix}}+\frac{E_{cut}}{{\sum }_{j\in B}{w}_{jx}}$$
A、B表示两个割集，并对所有节点的权重相加求和，对于有相同连接数的图像，是划分大小的一种粗糙度量方式。

D是对W每行元素求和过构成的对角矩阵，$d_i={\sum }_j{w}_{ij}$

归一化分割可通过最小化下面的优化问题而求得：
$$mi{n}_y\frac{y^T(D-W)y}{y^{T}Dy}$$

y包含离散标记，y只可以取特定的两个值，$y^{T}D$求和为0。

通过松弛约束条件让y取任意实数，最小化问题可以变为容易求解的特征分解问题，缺点是你需要对输出设定阈值或进行聚类，使它重新成为一个离散分割。

该问题成为：
$$L=D^{-1/2}W{D}^{-1/2}$$
难点在于定义边的权重。利用原始归一化后边的权重：

第二部分度量坐标矢量的接近程度。

from PIL import Image
from pylab import *
from numpy import *
from scipy.cluster.vq import *


def cluster(S,k,ndim):
    """ Spectral clustering from a similarity matrix."""
    
    # check for symmetry
    if sum(abs(S-S.T)) > 1e-10:
        print 'not symmetric'
    
    # create Laplacian matrix
    rowsum = sum(abs(S),axis=0)
    D = diag(1 / sqrt(rowsum + 1e-6))
    L = dot(D,dot(S,D))
    
    # compute eigenvectors of L
    U,sigma,V = linalg.svd(L,full_matrices=False)
    
    # create feature vector from ndim first eigenvectors
    # by stacking eigenvectors as columns
    features = array(V[:ndim]).T

    # k-means
    features = whiten(features)
    centroids,distortion = kmeans(features,k)
    code,distance = vq(features,centroids)
        
    return code,V


def ncut_graph_matrix(im,sigma_d=1e2,sigma_g=1e-2):
    """ Create matrix for normalized cut. The parameters are 
        the weights for pixel distance and pixel similarity. """
    
    m,n = im.shape[:2] 
    N = m*n
    
    # normalize and create feature vector of RGB or grayscale
    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 coordinates for distance computation
    xx,yy = meshgrid(range(n),range(m))
    x,y = xx.flatten(),yy.flatten()
    
    # create matrix with edge weights
    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归一化矩阵w中填充值。

第二个函数将拉普拉斯特征分解后的前ndim个特征向量合并在一起构成矩阵w，并对像素进行K-means聚类。

在样本图像上进行测试：

import ncut
from PIL import Image
from numpy import *
from pylab import *
im = array(Image.open('test/C-uniform33.ppm'))
m,n = im.shape[:2]
# 调整图像的尺寸大小为 (wid,wid)
wid = 50
rim = array(Image.fromarray(im))
rim.resize((wid,wid))
rim = 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)
code= array(Image.fromarray(code))
code.resize((m,n))

# 绘制分割结果
figure()
imshow(code)
gray()
show()

结果图：（可能因为版本的问题，结果出现分割的问题）

相似矩阵前4个特征向量

9.3 变分法

Chan-Vese 分割模型对于待分割图像区域假定一个分片常数图像模型。这里我们集中关注两个区域的情形，比如前景和背景，不过这个模型也可以拓展到多区域，

用一组曲线将图像分成两个区域，分割通过最小化模型能量给出：


若用$\lambda |c_1-c_2|$替换ROF(降噪)中的$\lambda$，则于ROF方程形式一致。

最小化Chan-Vese模型转变成设定阈值的ROF问题，调低阈值以确保足够的迭代次数：

import rof
from numpy import *
from PIL import Image

im = array(Image.open('houses.png').convert("L"))
U,T = rof.denoise(im,im,tolerance=0.001)
t = 0.4 # 阈值
from matplotlib.pyplot import *
imsave('result.pdf',U < t*U.max())
imshow(U < t*U.max(),'gray')