Python计算机视觉编程——第9章图像分割_python graph edge label-CSDN博客

本文链接：https://blog.csdn.net/fjyalzl/article/details/126504143

9.1 图割（Graph Cut）

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

9.1 图割（Graph Cut）

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

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

图割C（C是图中所有边的集合）的代价函数定义为所有割边的权重求和相加：

$\mathrm{E}_{\mathrm{cut}}=\sum_{(\mathrm{i}, \mathrm{j}) \in \mathrm{C}} \mathrm{w}_{\mathrm{ij}}$

wij是图中节点i到节点j的边（i，j）的权重，并且是对割C所以的边进行求和。

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

寻找最小割（minimum cut）等同于在源点和汇点间寻找最大流（maximum flow）。此外，很多高效的算法都可以解决这些最大流、最小割的问题。

我们在图割例子里将采用 python-graph 工具包，下载地址：https://github.com/pmatiello/python-graph 下载后将压缩文件解压到项目目录。代码运行时需要修改一些导入路径，否则可能会报错。

from pygraph.core.pygraph.classes.digraph import digraph
from pygraph.core.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）指一个像素与其正上方、正下方、左边、右边的像素直接相连。

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

下面给出创建这样一个图的步骤：

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

为边的权重建立如下模型：

$\begin{aligned} w_{s i} &=\frac{p_{F}\left(I_{i}\right)}{p_{F}\left(I_{i}\right)+p_{B}\left(I_{i}\right)} \\ w_{i t} &=\frac{p_{B}\left(I_{i}\right)}{p_{F}\left(I_{i}\right)+p_{B}\left(I_{i}\right)} \\ w_{i j} &=\kappa \mathrm{e}^{-\left|I_{i}-I_{j}\right|^{\prime} / \sigma} \end{aligned}$

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

创建一个名为 graphcut.py 的文件，将下面从一幅图像创建图的函数写入该文件中：

from pylab import *
from numpy import *

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

from PCV.classifiers import bayes

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


def build_bayes_graph(im, labels, sigma=1e2, kappa=1):
    """ 从像素四邻域建立一个图，前景和背景（前景用 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]) + 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:  # 左边存在
            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 x m + 2 的图。注意源点和汇点的索引；为了简化像素的索引，将最后的两个索引作为源点和汇点的索引。

为了在图像上可视化覆盖的标记区域，可以利用 contourf() 函数填充图像（这里指带标记图像）等高线间的区域，参数 alpha 用于设置透明度。将下面函数增加到 graphcut.py 中：

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')

图建立起来后便需要在最优位置对图进行分割。下面的函数可以计算最小割并将输出结果重新格式化为一个带像素标记的二值图像：

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 list(cuts.items())[:-2]:  # 不要添加源点 / 汇点
        res[pos] = label

    return res.reshape((m, n))

下面的例子会读取一幅图像，从图像的两个矩形区域估算出类概率，然后创建一个图：

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

im = array(Image.open('14.jpg'))

im = cv.resize(im, None, fx = 0.01, fy = 0.01, 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()

resize（）函数中数值决定了图像的缩放尺寸。根据选择图像不同，过大的尺寸可能会不适合 Python graph 库，导致报错。变量Kappa决定了近邻像素间的相对权重，随着K值增大，边界将变得更平滑，细节部分也逐渐丢失。

9.1.2 用户交互式分割

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

下面是一个例子，它会载入一幅图像及对应的标注信息，然后将其传递到我们的图割分割路径中：

# from scipy.misc import imresize
# from PCV.tools import graphcut
import graphcut
from PIL import Image
from pylab import *
# import numpy as np
# # import pylab as pl
# import matplotlib.pyplot as plt
def create_msr_labels(m, lasso=False):
    """ Create label matrix for training from
    user annotations. """
    # rim = im.reshape((-1,2))
    # m = m.convert("L")
    size = m.shape[:2]
    # m = Image.fromarray(m.astype('uint8')).convert("L")
    # size = m.shape[:2]
    labels = zeros(size)
    # background
    labels[m == 0] = -1
    labels[m == 64] = -1
    # foreground
    if lasso:
        labels[m == 255] = 1
    else:
        labels[m == 128] = 1
    return labels


# load image and annotation map
im = array(Image.open('1.png'))
m = array(Image.open('1-1.png').convert('L'))

# resize
# scale = 0.32
scale = 0.05  # scale = 0.32 ，scale*scale ~= 0.1跑起来非常慢,scale=0.05代码跑通比较快
# im = imresize(im, scale, interp='bilinear')
# m = imresize(m, scale, interp='nearest')
h1, w1 = im.shape[:2]
h2, w2 = m.shape[:2]
print(h1, w1)
print(h2, w2)
# num_px = int(h * np.sqrt(0.07))
# num_py = int(w * np.sqrt(0.07))
px1 = int(w1 * scale)
py1 = int(h1 * scale)
px2 = int(w2 * scale)
py2 = int(h2 * scale)
# imresize(im, 0.07,interp='bilinear')  ##imresize被scipy.misc弃用，用PIL库中的resize替代
im = array(Image.fromarray(im).resize((px1, py1), Image.BILINEAR))
m = array(Image.fromarray(m).resize((px2, py2), Image.NEAREST))
oim = im
print(im.shape[:2])
print(m.shape[:2])
# create training labels
labels = create_msr_labels(m, False)
print('labels finish')
# build graph using annotations
g = graphcut.build_bayes_graph(im, labels, kappa=2)
print('build_bayes_graph finish')
# cut graph
res = graphcut.cut_graph(g, im.shape[:2])
print('cut_graph finish')
# remove parts in background
res[m == 0] = 1
res[m == 64] = 1
# labels[m == 0] = 1
# labels[m == 64] = 1

# plot original image
fig = figure()
subplot(121)
imshow(im)
gray()
title('before')
axis('off')

# plot the result
subplot(122)
imshow(res)
gray()
xticks([])
yticks([])
title('after')
axis('off')

show()
fig.savefig('912.pdf')

print('finish')

实验发现，该算法对于画面构成简单地图像分割效果较好，且运行结果较快。当图像前景与背景相似（较为复杂）时，分割效果会显著下降，且运行时间较长。

9.2 利用聚类进行分割

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

$\mathrm{E}_{\mathrm{ncut}}=\frac{\mathrm{E}_{\mathrm{cut}}}{\sum_{\mathrm{i} \in \mathrm{A}} \mathrm{w}_{\mathrm{ix}}}+\frac{\mathrm{E}_{\mathrm{cut}}}{\sum_{\mathrm{j} \in \mathrm{B}} \mathrm{w}_{\mathrm{jx}}}$

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

定义 W 为边的权重矩阵，矩阵中的元素wij为连接像素 i 和像素 j 边的权重。D 为对 W 每行元素求和后构成的对角矩阵，即 $D=\operatorname{diag}\left(d_{i}\right)$ , $d_{i}=\sum_{j} W_{i j}$ 。归一化分割可以通过最小化下面的优化问题而求得:

$\min _{y} \frac{y^{T}(D-W) y}{y^{T} D y}$

向量y yy 包含的是离散标记，这些离散标记满足对于b 为常数 $y_{i} \in\{1,-b\}$ （即 y 只可以取这两个值）的约束， $y^{T}D$ 求和为 0。由于这些约束条件，该问题不太容易求解。

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

松弛该问题后，该问题便成为求解拉普拉斯矩阵特征向量问题：

$L=D^{-1 / 2} W D^{-1 / 2}$

正如谱聚类情形一样，现在的难点是如何定义像素间边的权重wij。我们利用原始归一化割论文中的边的权重，通过下式给出连接像素 i 和像素 j 的边的权重：

$w_{i j}=e^{-\left|I_{i}-I_{j}\right|^{2} / \sigma g} e^{-\left|x_{i}-x_{j}\right|^{2} / \sigma d}$

将下面的函数添加到名为 ncut.py 的文件中：

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

下面的脚本为一个完整的例子：

# from PCV.tools import ncut
# from scipy.misc import imresize
import ncut
from pylab import *
from PIL import Image

im = array(Image.open('uniform/train/C-uniform03.ppm'))
m, n = im.shape[:2]
print(n, m)
# resize image to (wid,wid)
wid = 50
# rim = imresize(im, (wid, wid), interp='bilinear')
rim = np.array(Image.fromarray(im).resize((wid, wid), Image.Resampling.BILINEAR))
rim = array(rim, 'f')
# create normalized cut matrix
A = ncut.ncut_graph_matrix(rim, sigma_d=1, sigma_g=1e-2)
# cluster
code, V = ncut.cluster(A, k=3, ndim=3)
print(array(V).shape)
print("ncut finish")

# 变换到原来的图像大小
# codeim = imresize(code.reshape(wid,wid),(m,n),interp='nearest')
codeim = array(Image.fromarray(code.reshape(wid, wid)).resize((n, m), Image.Resampling.NEAREST))
# imshow(imresize(V[i].reshape(wid,wid),(m,n),interp=’bilinear’))
# v = zeros((m,n,4),int)
v = zeros((4, m, n), int)
for i in range(4):
    v[i] = array(Image.fromarray(V[i].reshape(wid, wid)).resize((n, m), Image.Resampling.BILINEAR))

# 绘制分割结果
fig = figure()
gray()
subplot(242)
axis('off')
imshow(im)

subplot(243)
axis('off')
imshow(codeim)

for i in range(4):
    subplot(2, 4, i + 5)
    axis('off')
    imshow(v[i])

show()

因为 Numpy 的 linanlg.svd() 函数在处理大型矩阵时的计算速度并不够快（有时对于太大的矩阵甚至会给出不准确的结果），所以这里我们重新设定图像为一固定尺寸（在该例中为 50×50），以便更快地计算特征向量。在重新设定图像大小的时候我们采用了双线性插值法；因为不想插入类标记，所以在重新调整分割结果标记图像的尺寸时我们采用近邻插值法。注意，重新调整到原图像尺寸大小后第一次利用 reshape 方法将一维矩阵变换为（wid，wid）二维数组。

及时对于一些简单的例子，对图像进行阈值处理不会给出相同的结果，对RGB值或灰度值进行聚类也一样；原因是它们都没有考虑像素的近邻。

9.3 变分法

诸如 ROF 降噪、K-means 和 SVM 的例子，这些都是优化问题。当优化的对象是函数时，该问题称为变分问题，解决这类问题的算法称为变分法。我们看一个简单而有效的变分模型。

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

如果我们用一组曲线 $\Gamma$ 将图像分离成两个区域Ω1 和 Ω2，分割是通过最小化 Chan-Vese 模型能量函数给出的

$E(\Gamma)=\lambda \operatorname{length}(\Gamma)+\int_{\Omega 1}\left(I-c_{1}\right)^{2} d x+\int_{\Omega 2}\left(I-c_{2}\right)^{2} d x$

该能量函数用来度量与内部平均灰度常数 c1 和外部平均灰度常数 c2 的偏差。这里这两个积分是对各自区域的积分，分离曲线 $\Gamma$ 长度用以选择更平滑的方案。

由分片常数图像 $U=\chi_{1} c_{1}+\chi_{2} c_{2}$ ，我们可以将上式重写为：

$E(\Gamma)=\lambda \frac{\left|c_{1}-c_{2}\right|}{2} \int|\nabla U| \mathrm{d} \mathbf{x}+\|I-U\|^{2}$

import rof
from pylab import *
from PIL import Image
import imageio
from skimage import *

im1 = array(Image.open('7.jpg').convert("L"))
im2 = array(Image.open('4.png').convert("L"))
U1, T1 = rof.denoise(im1, im1, tolerance=0.001)
U2, T2 = rof.denoise(im2, im2, tolerance=0.001)

t1 = 0.8  # flower32_t0 threshold
t2 = 0.4  # ceramic-houses_t0 threshold
seg_im1 = img_as_uint(U1 < t1 * U1.max())
seg_im2 = img_as_uint(U2 < t2 * U2.max())

fig = figure()
gray()
subplot(231)
axis('off')
imshow(im1)

subplot(232)
axis('off')
imshow(U1)

subplot(233)
axis('off')
imshow(seg_im1)

subplot(234)
axis('off')
imshow(im2)

subplot(235)
axis('off')
imshow(U2)

subplot(236)
axis('off')
imshow(seg_im2)

show()

# scipy.misc.imsave('../images/ch09/flower32_t0_result.pdf', seg_im)
#imageio.imsave('../images/ch09/flower32_t0_result.pdf', seg_im1)
#imageio.imsave('../images/ch09/ceramic-houses_t0_result.pdf', seg_im2)
# fig.savefig('../images/ch09/flower32_t0_result.pdf', seg_im1)
# fig.savefig('../images/ch09/ceramic-houses_t0_result.pdf', seg_im2)