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.shapeunary = np.zeros([row,col,2])for x in range(0,col):for y in range(0,row):# pixel colorpc = im[y,x,:]# data termfdf = dis(pc,fg)fdb = dis(pc,bg)# prior term# right neighbor pixelfpr = 1# below neighbor pixelfpb = 1unary[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.shapesmooth = 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.shapedim = np.zeros(im.shape)for i in range(0,row):for j in range(0,col):# backgroundif labels[i,j] == 1:dim[i,j] = bgelif labels[i,j] == 0:dim[i,j] = fgreturn 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 gcoim1 = 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))/3def median(data):data.sort()half = len(data) // 2return (data[half] + data[~half])/2def mode(data):data = list(data)return max(set(data),key=data.count)

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

def Ddistribution(im1,im2):size1 = im1.shapesize2 = im2.shapeif 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) / 3disparity = np.argmin(dis)if dis[disparity] < 5 :Dis[j] = disparitydisp = [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.shapesize2 = im2.shapeif 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 colorpixel1 = 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 / 255unary[x,y] = labels#         print('label computing: %.2f%%'%(x/row)*100)unary_cut = unary[:,:int(col-dmax),:] * depthreturn unary_cut, depth

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

# graph cut
def gcut(unary,depth,lmd):[row,col,c] = unary.shapefp = 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 labelsDisparity = 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和数据项平均值的比值很重要。