【python计算机视觉编程——9.图像分割】

9.图像分割

9.1 图割

可以选择不装Graphviz，因为原本觉得是要用，后面发现好像用不到。不安装可直接跳到下一步

安装Graphviz

1. 首先需要先下载Graphviz软件（Download | Graphviz），那些包先不要下载，网上说先下载包再下载软件会报错。在安装过程中，需要注意下图中的一步，其余都是一直下一步就行
   

2. 检查一下环境变量的路径
   

3. 接着在自己创建的虚拟环境下安装包
   pip install pydotplus

   pip install graphviz

4. 这里需要注意的是，还需要再安装一个包，否则单单安装上面的会报错
   pip install python-graphviz

5. 测试代码

   from graphviz import Digraph
   dot = Digraph(comment='The Round Table')
   dot.node('A', 'King Arthur')
   dot.node('B', 'Sir Bedevere the Wise')
   dot.node('L', 'Sir Lancelot the Brave')
   
   dot.edges(['AB', 'AL'])
   dot.edge('B', 'L', constraint='false')
   print(dot.source)  
   
   dot.render('round-table.gv',format='jpg', view=True)  
   

   

   

下一步：正文

• 图割：将一个有向图分割成两个互不相交的集合
• 基本思想：相似且彼此相近的像素应该划分到同一区域

图割C（C是图中所有边的集合）的“代价”函数定义为所有割的边的权重求合相加：
$$E_{cut}=\sum _{(i,j)\in C}{w}_{ij}$$
$w_{ij}$是图中节点i到节点j的边$(i,j)$的权重，并且是对割C所有的边进行求和

我们需要用图来表示图像，并对图进行划分，以使得$E_{cut}$最小。同时在用图表示图像时，需要额外增加两个节点（源点和汇点），并仅考虑那些将源点和汇点分开的割

寻找最小割等同于在源点和汇点间寻找最大流，这里需要用到python-graph工具包（$注：不是pip下载!$），工具包地址如下：GitHub - pmatiello/python-graph: New official repository: https://github.com/Shoobx/python-graph

下载完后，把文件夹放入导包的根目录

根据路径进行引包，如果引入没有报错，就说明没有问题

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

这里我是报错了：“jaraco.text"中没有drop_comment, join_continuation, yield_lines函数的问题，然后我在”_jaraco_text.py"文件里找到了这三个函数，索性就直接把他提到根目录上，发现就没报错了

另一个导包路径错误在"digraph.py"和"minmax.py"文件中

接着就可以运行代码了

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)

9.1.1 从图像创建图

我们需要利用图像像素作为节点定义一个图，除了像素节点外，还有两个特定的节点——“源”点和“汇”点，来分别代表图像的前景和背景，我们需要做的是将所有像素与源点、汇点链接起来。

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

接着需要用朴素贝叶斯分类器进行分类，我们将第8章的BayesClassifier类搬过来

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 = 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] / (np.linalg.norm(vim[i]) + 1e-9)
    
    # go through all nodes and add edges
    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*np.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*np.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*np.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*np.exp(-1.0*sum((vim[i]-vim[i+n])**2)/sigma)
            gr.add_edge((i,i+n),wt=edge_wt)
    return gr    

def gauss(m,v,x):
    """ Evaluate Gaussian in d-dimensions with independent 
        mean m and variance v at the points in (the rows of) x. 
        http://en.wikipedia.org/wiki/Multivariate_normal_distribution """
    
    if len(x.shape)==1:
        n,d = 1,x.shape[0]
    else:
        n,d = x.shape
            
    # covariance matrix, subtract mean
    S = np.diag(1/v)
    x = x-m
    # product of probabilities
    y = np.exp(-0.5*np.diag(np.dot(x,np.dot(S,x.T))))
    
    # normalize and return
    return y * (2*np.pi)**(-d/2.0) / (np.sqrt(np.prod(v)) + 1e-6)

写入新函数

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 = BayesClassifier()
    bc.train(train_data)

    # 获取所有像素的概率
    bc_lables,prob = bc.classify(vim)
    prob_fg = prob[0]
    prob_bg = prob[1]
    
    # 用m*n+2 个节点创建图
    gr = nx.DiGraph()
    
    nodes=[]
    for i in range(m*n+2):
        nodes.append(str(i))
    gr.add_nodes_from(nodes)
    
    
    source = m*n #  倒数第二个是源点
    sink = m*n+1 # 最后一个节点是汇点

    #  归一化
    for i in range(vim.shape[0]):
        vim[i] = vim[i] / (np.linalg.norm(vim[i]) + 1e-9)
    
    # go through all nodes and add edges
    for i in range(m*n):
        # 从源点添加边
        gr.add_edge(str(source),str(i),capacity=prob_fg[i]/(prob_fg[i]+prob_bg[i]))
        
        # 向汇点添加边
        gr.add_edge(str(i),str(sink),capacity=prob_bg[i]/(prob_fg[i]+prob_bg[i]))
        
        # 向相邻节点添加边
        if i%n != 0: # 左边存在
            edge_wt = kappa*np.exp(-1.0*sum((vim[i]-vim[i-1])**2)/sigma)
            gr.add_edge(str(i),str(i-1),capacity=edge_wt)
        if (i+1)%n != 0: # 如果右边存在
            edge_wt = kappa*np.exp(-1.0*sum((vim[i]-vim[i+1])**2)/sigma)
            gr.add_edge(str(i),str(i+1),capacity=edge_wt)
        if i//n != 0: # 如果上方存在
            edge_wt = kappa*np.exp(-1.0*sum((vim[i]-vim[i-n])**2)/sigma)
            gr.add_edge(str(i),str(i-n),capacity=edge_wt)
        if i//n != m-1: # 如果下方存在
            edge_wt = kappa*np.exp(-1.0*sum((vim[i]-vim[i+n])**2)/sigma)
            gr.add_edge(str(i),str(i+n),capacity=edge_wt)
    return gr    

def show_labeling(im,labels):
    """显示图像的前景和背景区域。前景labels=1,背景labels=-1，其他labels＝0 """
    
    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):
    """    Solve max flow of graph gr and return binary 
        labels of the resulting segmentation."""
#     print(gr)
    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)
#     print(cuts)
    # convert graph to image with labels
    res = np.zeros(m*n)
    for pos,label in list(cuts.items())[:-2]: # 遍历所有节点，忽略源节点和汇节点
        # 但因为cuts.items()返回的是元组，需先转成列表再进行切片
        res[pos] = label

    return res.reshape((m,n))

其中书本中from scipy.misc import imresize模块，已经不存在于imresize中，这里使用Pillow库中的resize函数进行替代 resize_image_pillow

def resize_image_pillow(image_path, output_path, scale_factor):
    # 打开图像文件
    with Image.open(image_path) as img:
        # 计算新的尺寸
        new_width = int(img.width * scale_factor)
        new_height = int(img.height * scale_factor)
        
        # 使用双线性插值调整图像大小
        img_resized = img.resize((new_width, new_height), resample=Image.BILINEAR)
        
        # 保存调整后的图像
#         return img_resized
        img_resized.save(output_path)

import numpy as np
from PIL import Image
from pylab import *

# resize_image_pillow('empire.jpg', 'empire.jpg', 0.07)
im=np.array(Image.open('empire.jpg'))
size=im.shape[:2]
labels=np.zeros(size)
labels[3:18,3:18]=-1
labels[-18:-3,-18:-3]=1

# 对图进行分割
g = build_bayes_graph(im,labels,kappa=1)
res=cut_graph(g,size)

figure()
show_labeling(im,labels)

figure()
imshow(res)
gray()
axis('off')

show()

9.1.2 用户交互式分割

def create_msr_labels(m, lasso=False):
    """ Create label matrix for training from
    user annotations. """
    labels = np.zeros(im.shape[:2])
    # 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('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 = build_bayes_graph(im, labels, kappa=2)

# cut graph
res = 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')

9.2 利用聚类进行分割

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

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,full_matrices=False)

    # 从前 ndim 个特征向量创建特征向量
    # 堆叠特征向量作为矩阵的列
    features = array(V[:ndim]).T

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

在运行下面代码之前，需要安装scikit-image，记得在自己的虚拟环境下安装（我用pip安装不了，后面改用conda，只要在虚拟环境下，用哪个（pip或conda）都是安装在虚拟环境下）

conda install scikit-image

import cv2
import numpy as np
from pylab import  *
from PIL import Image
from skimage.transform import resize

im = Image.open('empire.jpg')
m,n = np.array(im).shape[:2]
# 调整图像的尺寸大小为(wid,wid)
wid = 50


rim = im.resize((50,50),Image.BILINEAR)
rim = array(rim,'f')
# 创建归一化割矩阵
# print(rim.shape[:2] )
A = ncut_graph_matrix(rim,sigma_d=1,sigma_g=1e-2)
# 聚类
code,V=cluster(A,k=3,ndim=3)
# 变换到原来的图像大小


image=code.reshape(wid,wid)
print(image)


codeim = resize(image,(m,n),mode='reflect',anti_aliasing=False,order=0)
# 绘制分割结果
figure()
imshow(codeim)
gray()
show()

9.3 变分法

当优化的对象是函数时，该问题称为变分问题，需要使用ROF进行降噪。

denoise函数需要传入以下参数

• im: 输入的噪声图像（灰度图像）。
• U_init: 对 U（去噪图像）的初始猜测。
• tolerance: 收敛的容忍度，用于判断迭代是否结束。
• tau: 步长（或称为步伐），用于控制更新的幅度。
• tv_weight: 总变差正则化项的权重，控制去噪程度。

denoise函数返回参数

• U: 去噪后的图像。
• im - U: 图像的纹理残差，即原始图像中未被去噪部分的残余。

def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):
    """ 这个函数实现了 Rudin-Osher-Fatemi (ROF) 去噪模型，
    ROF 模型是一个常用的图像去噪方法，基于总变差(Total Variation, TV)
    正则化来去除噪声，同时保留图像的边缘信息"""
        
    m,n=im.shape  #获取图像的高度和宽度
    
    #初始化
    U=U_init
    Px=im      # 对偶域的x分量
    Py=im      # 对偶域的y分量
    error=1
    while(error>tolerance):
        Uold=U

        #原始变量的梯度
        GradUx=roll(U,-1,axis=1)-U       #变量U梯度的x分量
        GradUy=roll(U,-1,axis=0)-U       #变量U梯度的y分量

        #更新对偶变量
        PxNew=Px+(tau/tv_weight)*GradUx  #更新Px
        PyNew=Py+(tau/tv_weight)*GradUy  #更新Py
        NormNew=maximum(1,sqrt(PxNew**2+PyNew**2))#计算PxNew和PyNew的范数,确保其最小值为1

        
        Px=PxNew/NormNew    #更新x分量
        Py=PyNew/NormNew    #更新y分量
        
        RxPx=roll(Px,1,axis=1)#计算Px在x方向上的右移
        RyPy=roll(Py,1,axis=0)#计算Px在y方向上的下移
        
        DivP=(Px-RxPx)+(Py-RyPy)#计算Px和Py的梯度
        
        U=im+tv_weight*DivP # 更新去噪后的图像U
        
        error=linalg.norm(U-Uold)/sqrt(n*m)# 计算当前误差
    return U,im-U            #返回去噪后的图像U和噪声图像

因为 scipy.misc.imsave 已被弃用，所以需要用其他库来完成，这里使用Pillow库来保存图像

import numpy as np
from PIL import Image
im = np.array(Image.open('ceramic-houses_t0.png').convert('L'))
U,T=denoise(im,im,tolerance=0.001)
t=0.4


# 基于阈值生成二值图像
binary_image = U < t * U.max()

# 将布尔数组转换为 uint8 格式（0 或 255）
binary_image_uint8 = (binary_image * 255).astype(np.uint8)

# 创建 Image 对象
img = Image.fromarray(binary_image_uint8)

# 保存图像为 PDF
img.save('result.pdf')

from pylab import  *
gray()
subplot(121)
imshow(U)
subplot(122)
imshow(img)