マルコフ確率場上で反復条件付きモードによるノイズ除去を組んでみる(Python)

反復条件付きモードとは

反復条件付きモードはICMともよばれ，マルコフ確率場上で各ノードの状態を順番に変更しながら，ノード全体の状態を最適化するためのアルゴリズムです．しかしこのアルゴリズムはノード毎に勾配法を適用するだけなので，大域解に収束せず局所解しか求められません．

使用データ

ソースコード解説

大まかな流れは，各ノードにおいて，そのノードがとりうる状態全てについてエネルギーを求めてあげます．このエネルギーはそのノードがどの状態をとるべきかの度合いを表し，各ノードはエネルギーが最も低くなる状態をとります．

まずreference変数にノイズ画像をコピーします．これがノードの状態になります．次に各ノードについて隣接するノードの状態を調べながら自分がどの状態をとるべきかのエネルギーを以下の様に計算します．

for i in range(1, height - 1):
    for j in range(1, width - 1):
        for k in range(len(dy)):

            #近傍のノードの状態を考慮する
            if reference[i+dy[k],j+dx[k]] == 1:
               prob1[i,j] -= beta
               prob0[i,j] += beta
            else:
               prob1[i,j] += beta
               prob0[i,j] -= beta

次に観測値(ノイズ画像)のみを考慮して各ノードがどの値をとるべきかのエネルギーを計算してあげます．それが以下になります．

prob1[image == 1] -= alpha
prob0[image == 0] -= alpha

ソースコード全体は以下の様になります．

#encoding:utf-8
import numpy as np
import cv2
import matplotlib.pyplot as plt
from skimage.filters import threshold_otsu

def MRF(image):
    iter = 30
    height, width = image.shape
    dy = [-1,-1,-1,0,0,1,1,1]
    dx = [-1,0,1,-1,1,-1,0,1]
    reference = np.copy(image)
    alpha = 2.5
    beta = 4.0

    for t in range(iter):
        print(t)

        prob0 = np.zeros(reference.shape)
        prob1 = np.zeros(reference.shape)

        for i in range(1, height - 1):
            for j in range(1, width - 1):
                for k in range(len(dy)):

                    #近傍のノードの状態を考慮する
                    if reference[i+dy[k],j+dx[k]] == 1:
                        prob1[i,j] -= beta
                        prob0[i,j] += beta
                    else:
                        prob1[i,j] += beta
                        prob0[i,j] -= beta

        #観測値を考慮する
        prob1[image == 1] -= alpha
        prob0[image == 0] -= alpha

        reference[prob1 > prob0] = 0
        reference[prob1 < prob0] = 1

    return reference


def addNoise(image):
    output = np.copy(image)
    flags  = np.random.binomial(n=1, p=0.05, size=image.shape)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            if flags[i,j]:
                output[i,j] = not(output[i,j])

    return output

def main():
    image = cv2.imread("Lenna.png", 0)
    binary = image > threshold_otsu(image)
    noise = addNoise(binary.astype(np.int))
    denoise = MRF(noise)

    plt.gray()
    plt.subplot(131)
    plt.imshow(image)
    plt.subplot(132)
    plt.imshow(noise)
    plt.subplot(133)
    plt.imshow(denoise)
    plt.show()

if __name__ == '__main__':
    main()

実行結果

確かにノイズが除去されています．ＩＣＭは非常に簡単に組めるのでマルコフ確率場の入門にはお勧めです．