2024.9.14 Python与图像处理新国大EE5731课程大作业，马尔可夫随机场和二值图割，校正立体图像的深度

1.马尔科夫随机场和二值图割

马尔可夫随机场（MRF, Markov Random Field）：
MRF 是一种用来描述图像像素之间空间关系的概率模型。它假设图像中的像素不仅取决于自身的值，还与周围像素有关。这种模型经常用于图像分割、去噪等任务。
在去噪问题中，MRF 可以用于表示像素之间的关联性，确保去噪过程中不仅关注单个像素，还考虑周围像素的影响。
二值图割（Binary Graphcuts）：
图割算法是一种常用于图像分割和去噪的技术。它通过将图像像素建模为图中的节点，并使用图割技术来最小化能量函数，从而实现分割或去噪。
二值图割的意思是将像素分类为两个类别，通常是“前景”和“背景”，或者“噪声”和“非噪声”。
结合图论中的最小割问题，可以找到一种最优的像素分割方式，从而去除噪声。

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
import random as rm
import math
import cv2
import random
import gco

# foreground blue
fg = np.array([0,0,255])
# background yellow
bg = np.array([245,210,110])

im = cv2.imread('bayes_theorem.jpg',cv2.IMREAD_COLOR)
im = cv2.cvtColor(im, cv2.COLOR_BGR2RGB)

#show the original figure
plt.figure( figsize=(15,15) )
plt.imshow(im)

# distance
def dis(xi,xj):
    return np.sum(np.abs(xi - xj))/3/255

# generate nodes
def nodes(im,lmd,fg,bg):
    
    [row,col,c] = im.shape
    unary = np.zeros([row,col,2])
    for x in range(0,col):
        for y in range(0,row):

            # pixel color
            pc = im[y,x,:]

            # data term
            fdf = dis(pc,fg)
            fdb = dis(pc,bg)

            # prior term
            # right neighbor pixel
            fpr = 1
            # below neighbor pixel
            fpb = 1

            unary[y,x,0] = fdf + lmd*(fpr + fpb)
            unary[y,x,1] = fdb + lmd*(fpr + fpb)
    
    return unary

# graph cut
def gcut(unary,lmd):
    [row,col,c] = unary.shape
    smooth = 1 - np.eye(2)
    labels = gco.cut_grid_graph_simple(unary, smooth*lmd, n_iter=-1)
    labels = labels.reshape(row,col)
    return labels

# original iamge denoise
def dimage(im,labels):
    [row,col,c] = im.shape
    dim = np.zeros(im.shape)
    for i in range(0,row):
        for j in range(0,col):

            # background
            if labels[i,j] == 1:
                dim[i,j] = bg
            elif labels[i,j] == 0:
                dim[i,j] = fg

    return dim.astype(int)

# lambda = 1
lmd = 1
# nodes
unary = nodes(im,lmd,fg,bg)
# graphcut
label = gcut(unary,lmd)
# image denoising
dim = dimage(im,label)

这个代码实现了一个使用图割算法（Graph Cut）进行图像去噪的简单例子。代码的主要思想是使用马尔可夫随机场（MRF）模型，将图像像素分为两类：前景（蓝色）和背景（黄色），并通过图割来优化像素的分割，以达到去除噪声的目的。
1.距离函数dis():这个函数用于计算两个像素点之间的色彩距离。这里计算的是每个像素的绝对色差，然后进行归一化处理（将颜色差值除以 255，保持在 0-1 之间）
2.生成节点：nodes()，这个函数的作用是构建每个像素的代价（unary term），即每个像素分别属于前景和背景的代价。
3.gcut()图割算法通过最小化能量函数，分割前景和背景，返回每个像素的标签（0：前景，1：背景）在这一步中已经进行了最小化处理了，给每个像素一个前景或者后景的标签
4.def dimage(im,labels):这个函数染色
5.通过设置 lambda 参数来调整平滑项的权重。较大的 lambda 会增强对像素间平滑性的惩罚，使图像看起来更平滑，但也可能会损失细节。生成 unary 代价矩阵，执行图割，最终生成去噪后的图像。

效果图：

2.校正立体图像的深度

主要任务是图像匹配和视差计算，用来将两个深度图像进行对比，从而估计深度信息。具体来说，通过计算两幅图像中对应像素点之间的匹配，生成一个视差图，进而可能用于3D重建或者深度信息提取。

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
import random as rm
import math
import cv2
import random
import gco

im1 = cv2.imread('depth_im1.png',cv2.IMREAD_COLOR)
im1 = cv2.cvtColor(im1, cv2.COLOR_BGR2RGB)
im2 = cv2.imread('depth_im2.png',cv2.IMREAD_COLOR)
im2 = cv2.cvtColor(im2, cv2.COLOR_BGR2RGB)

im = np.hstack((im1,im2))
#show the original figure
plt.figure( figsize=(15,15) )
plt.imshow(im)

使用 OpenCV 读取两张深度图（im1 和 im2），并将它们水平拼接成一张图，然后显示出来。

# distance
def distance(xi,xj):
    return np.sum(np.abs(xi - xj))/3

def median(data):
    data.sort()
    half = len(data) // 2
    return (data[half] + data[~half])/2

def mode(data):
    data = list(data)
    return max(set(data),key=data.count)

distance：计算两个像素的欧几里得距离
median：返回列表的中值
mode：返回列表的众数值

def Ddistribution(im1,im2):
    
    size1 = im1.shape
    size2 = im2.shape
    
    if size1 != size2 : 
        raise ValueError('input shape not matches')
        
    dis_col = []
    for i in range(size1[0]):
        Dis = np.zeros(size1[1])
        for j in range(size1[1]):
            pixel = im1[i,j]
            im2_pixel = im2[i,j:]
            difference = abs(im2_pixel - pixel)
            dis = np.sum(difference,axis = 1) / 3
            disparity = np.argmin(dis)
            if dis[disparity] < 5 :
                Dis[j] = disparity
        
        disp = [p for p in Dis if p > 0]
        dis_col.append(mode(disp))
                        
    return dis_col

分布：遍历一幅图像中的所有像素，得到两幅图像的视差值范围
计算图像 im1 和 im2 之间的视差分布，逐行计算每个像素在 im1 和 im2 之间的最小差异，并且记录最匹配的位移量（视差）。

# generate nodes
def nodes(im1,im2,Disparity,depth = None):
    
    # input protection 
    size1 = im1.shape
    size2 = im2.shape
    
    if size1 != size2 : 
        raise ValueError('input shape not matches')
        
    [row,col,c] = size1
    # preprocessing: get the disparity distribution first for getting D: dmin and dmax 
#     Disparity = Ddistribution(im1,im2)
    dmin = min(Disparity) 
    dmax = max(Disparity) 
    if depth == None:
        depth = int(dmax - dmin)
    else :
        depth = int(depth)
    step = (dmax - dmin) / depth
#     print(Disparity)
    print('dmin:',dmin,'dmax',dmax)
    unary = np.zeros([row,col,depth])
    
    for x in range(0,row):
        for y in range(0,int(col - dmax)):

            # pixel color
            pixel1 = im1[x,y]
            labels = np.zeros(depth)
            
            for i in range(depth):
                
                if i < depth - 1:
                    pixel2_index_start = round(i * step + dmin)
                    pixel2_index_end = round((i+1) * step + dmin)
                elif i == depth - 1:
                    pixel2_index_start = round(i * step + dmin)
                    pixel2_index_end = int(dmax)
                
                pixel2 = im2[x,y+pixel2_index_start:y+pixel2_index_end]
                dis = []
                for p in pixel2:
                    dis.append(distance(pixel1,p))
                dis_min = min(dis)
                
                labels[i] = dis_min / 255
            
            unary[x,y] = labels
        
#         print('label computing: %.2f%%'%(x/row)*100)
    
    unary_cut = unary[:,:int(col-dmax),:] * depth
    return unary_cut, depth

节点：根据计算出的视差生成数据项，D(dmax - dmin)的深度等于视差的最大值减去视差的最小值
生成图像匹配的能量（unary）矩阵，每个像素点有多个视差候选值（基于不同的深度），该函数计算每个像素与可能的匹配点之间的差异并归一化为能量值。
Disparity 用于确定视差的上下限。
作用: 这个步骤通过遍历图像中的每一个像素点，计算它在不同视差下的匹配质量，输出的 unary 是图像的视差匹配成本图，后续会用于图割算法。

# graph cut
def gcut(unary,depth,lmd):
    [row,col,c] = unary.shape
    fp = np.zeros([depth,depth])
    for i in range(depth):
        for j in range(depth):
                fp[i,j] = abs(i-j)
    labels = gco.cut_grid_graph_simple(unary, fp*lmd, connect = 8, n_iter=-1)
    labels = labels.reshape(row,col)
    labels = labels[:,depth:]
    return labels

Disparity = Ddistribution(im2,im1)
print(min(Disparity),max(Disparity))

unary,depth = nodes(im2,im1,Disparity)

gcut：先生成先验项，然后进行图切割并返回标签
使用图割算法（graph cut）对能量矩阵进行全局优化，最终得到每个像素点的最佳视差值（标签）。这里的图割通过计算视差标签之间的平滑约束，确保邻近的像素视差值不会有太大跳变。
作用: 在视差图中，优化每个像素点的视差标签，使得输出的结果更连贯和自然。

lmd = 0.1
labels = gcut(unary,depth,lmd)

plt.figure( figsize=(8,8) )
plt.imshow(labels,'gray')

主要逻辑和输出
读取图像：首先读取并拼接两张图像，用于可视化差异。
计算视差分布：通过 Ddistribution 函数计算 im1 和 im2 之间的像素差异，提取视差分布。
生成能量矩阵：通过 nodes 函数计算能量矩阵，每个像素点在不同视差下的匹配成本。
图割优化：通过 gcut 函数使用图割算法对能量矩阵进行全局优化，得到最佳视差标签。
可视化视差图：最终输出视差标签，并可视化为灰度图。
lmd=0.1噪声大，锯齿明显

lmd=0.5鲁棒性更好

lmd=5

现在手动设定depth为10，之前的depth差不多是50左右

深度图像不明显。
depth 代表的是离散的深度层数，用于分离图像中不同像素的深度信息。它决定了在图像的视差计算过程中，离散化深度的精细程度。depth 的计算方式与视差范围有关，即两幅图像之间像素点位移的可能值范围（最小视差到最大视差之间的距离）

现在通过更改lmd的值可以看到不同的效果，lmd小的时候，噪声很大，lmd大的时候，深度信息就看不清了，对比相同lambda，不同dpeth的深度图，在dmin和dmax不变的情况下，深度层数越少，深度图中的噪声越小。同时lambda值和数据项、先验项的值有关，其实lambda和数据项平均值的比值很重要。