獨立隨機序列和為馬爾可夫序列，馬爾可夫隨機場的python實現

? ? ? ?前幾天看圖網絡的時候，不知怎么從貝葉斯網絡突然跳到了馬爾科夫隨機場，感覺還有點意思，不過還是沒有完全理清其中的邏輯，網上講的比較亂,參考博主on2way這篇博文《從貝葉斯理論到圖像馬爾科夫隨機場》我梳理了一下思路，找時間總結一下，先挖一個坑吧，由于博主給出的是matlab實現，我本身不怎么使用matlab，于是就簡單的用python復現了一下思路，由于沒有使用kmeans算法初始化，結果差距比較明顯，但是還算不錯，其中應該是有一些問題，包括最后多分類的接口沒有調好，還有四分類顯示結果貌似不太對，由于后天要考英語今天看來是沒法調完了，依舊是放著這個坑，有時間的時候代碼改一下，以下是python3的代碼，直接用jupyter notebook寫的：

import numpy as np
import cv2
import copy
from scipy.stats import norm
import matplotlib as plotimg = cv2.imread('Lena.jpeg',0)
cluster_num = 4
maxiter = 50
#隨機初始化標簽
label = np.random.randint(0,cluster_num,img.shape)def cal_pp(label):result = np.zeros((cluster_num, 8, label.shape[0], label.shape[1]))statis_result = np.zeros((cluster_num, label.shape[0], label.shape[1]))for i in range(cluster_num):label_temp = copy.deepcopy(label)#先設置為4再設置為0，防止處理0時出錯label_temp[label_temp!=i] = 4label_temp[label_temp==i] = 1label_temp[label_temp==4] = 0#根據label_temp可以計算出8個方向的標簽#從左上角起順時針方向，依次為0-8result[i][0][:,0] = 0result[i][0][0,:] = 0result[i][0][1:,1:] = label_temp[:-1,:-1]result[i][1][:,0] = 0result[i][1][1:,:] = label_temp[:-1,:]result[i][2][:,0] = 0result[i][2][0,:] = 0result[i][2][1:,1:] = label_temp[1:,:-1]result[i][3][:,0] = 0result[i][3][:,1:] = label_temp[:,:-1]result[i][4][:,-1] = 0result[i][4][:,:-1] = label_temp[:,1:]result[i][5][:,0] = 0result[i][5][-1,:] = 0result[i][5][1:,:-1] = label_temp[1:,:-1]result[i][6][-1,:] = 0result[i][6][:-1,:] = label_temp[1:,:]result[i][7][:,-1] = 0result[i][7][-1,:] = 0result[i][7][:-1,:-1] = label_temp[1:,1:]statis_result[i] = result[i][0]+result[i][1]+result[i][2]+result[i][3]+result[i][4]+result[i][5]+result[i][6]+result[i][7]statis_result = np.array(statis_result,dtype=np.float) / 8#防止出現0，原因不清楚statis_result[statis_result==0] = 0.001return statis_resultdef cal_lf(img, label):result = np.zeros((cluster_num, 8, label.shape[0], label.shape[1]))distribution_parameter = np.zeros((cluster_num,2),np.float)for i in range(cluster_num):img_temp = copy.deepcopy(img)img_temp[label!=i] = 0#根據原圖灰度統計每一個標簽的分布和方差img_list = img_temp.tolist()img_list = [item for sublist in img_list for item in sublist if item!= 0]#再根據原圖的灰度數值算出似然函數值distribution_parameter[i][0] = np.mean(img_list)distribution_parameter[i][1] = np.std(img_list)#計算每一點屬于不同類型的概率result_lf = np.zeros((cluster_num,label.shape[0],label.shape[1]),np.float)for i in range(cluster_num):result_lf[i] = norm(distribution_parameter[i][0],distribution_parameter[i][1]).pdf(img)return result_lfdef update_label(img, init_label, maxiter):label = init_labelfor i in range(maxiter):prior_probability = cal_pp(label)likelihood = cal_lf(img,label)probability = np.zeros((cluster_num,label.shape[0],label.shape[1]))for j in range(cluster_num):probability[j] = prior_probability[j]*likelihood[j]#根據概率值更新labelmax_probability = probability.max(0)#label_one不需要，這一部分代碼需要做一下調整，時間緊迫先演示一下label_two = probability[1]-max_probabilitylabel_two[label_two!=0] = 4label_two[label_two==0] = 1label_two[label_two==4] = 0label_three = probability[2]-max_probabilitylabel_three[label_three!=0] = 4label_three[label_three==0] = 2label_three[label_three==4] = 0label_four = probability[3]-max_probabilitylabel_four[label_four!=0] = 4label_four[label_four==0] = 3label_four[label_four==4] = 0label= label_two + label_three + label_four#每次繪制一張圖片cv2.imshow("test",label)cv2.waitKey(500)cv2.destroyAllWindows()return label#main函數中運行，博主使用jupyter notebook就沒有改了，請需要的自行改一下
final = update_label(img,label,maxiter)

原圖放這里供大家下載使用：

貼一下大致效果：
其中前兩排是第1-8次迭代結果，最后一排是第12、16、20、24次迭代結果，效果還行