# 用 Python 通过马尔可夫随机场（MRF）与 Ising Model 进行二值图降噪

### 建模

• yi：噪声图中的像素
• xi：原图中的像素，对应噪声图中的 yi

yi和xi有很强的联系。

xi和与它相邻的其他像素也存在较强的联系

1. {xi, yi}，即原图像素与噪声图像素对
2. {xi, xj}，其中 xj 表示与 xi 相邻的像素

ψC(xC) 即所谓的 potential function，为了方便我们通常只求它的对数形式 E(xC)（按照其物理意义称为 energy function

import numpy as np

def E_generator(beta, eta, h):
"""Generate energy function E.

Usage: E = E_generator(beta, eta, h)
Formula:
E = h * \sum{x_i} - beta * \sum{x_i x_j} - eta * \sum{x_i y_i}
"""
def E(x, y):
"""Calculate energy for matrices x, y.

Note: the computation is not localized, so this is quite expensive.
"""
# sum of products of neighboring paris {xi, yi}
xxm = np.zeros_like(x)
xxm[:-1, :] = x[1:, :]  # down
xxm[1:, :] += x[:-1, :]  # up
xxm[:, :-1] += x[:, 1:]  # right
xxm[:, 1:] += x[:, :-1]  # left
xx = np.sum(xxm * x)
xy = np.sum(x * y)
xsum = np.sum(x)
return h * xsum - beta * xx - eta * xy

return E

# sum of products of neighboring paris {xi, yi}
xxm = np.empty_like(x)
xxm[:-1, :] = x[1:, :]  # down
xxm[1:, :] += x[:-1, :]  # up
xxm[:, :-1] += x[:, 1:]  # right
xxm[:, 1:] += x[:, :-1]  # left
xx = np.sum(xxm * x)

### Iterated Conditional Modes

def E_generator(beta, eta, h):

def E(x, y):
...  # as shown before

def is_valid(i, j, shape):
"""Check if coordinate i, j is valid in shape."""
return i >= 0 and j >= 0 and i < shape[0] and j < shape[1]

def localized_E(E1, i, j, x, y):
"""Localized version of Energy function E.

Usage: old_x_ij, new_x_ij, E1, E2 = localized_E(Ecur, i, j, x, y)
"""
oldval = x[i, j]
newval = oldval * -1  # flip
# local computations
E2 = E1 - (h * oldval) + (h * newval)
E2 = E2 + (eta * y[i, j] * oldval) - (eta * y[i, j] * newval)
adjacent = [(0, 1), (0, -1), (1, 0), (-1, 0)]
neighbors = [x[i + di, j + dj] for di, dj in adjacent
if is_valid(i + di, j + dj, x.shape)]
E2 = E2 + beta * sum(a * oldval for a in neighbors)
E2 = E2 - beta * sum(a * newval for a in neighbors)
return oldval, newval, E1, E2

return E, localized_E

localized_Ex[i, j] 更改为另一状态，减去原来状态在 E 里对应的那几项，计算新状态对应的项，再加上去。adjacentneighbors均是为了过滤掉非法的边界元素而设。

ICM 的实现如下，为了方便画图查看 energy 变化情况，在每遍历完一行的时候，我们可以记录当前的时间和 energy：

def ICM(y, E, localized_E):
"""Greedy version of simulated_annealing()."""
x = np.array(y)
Ebest = Ecur = E(x, y)  # initial energy
initial_time = time.time()
energy_record = [[0.0, ], [Ebest, ]]

accept, reject = 0, 0
for idx in np.ndindex(y.shape):  # for each pixel in the matrix
old, new, E1, E2 = localized_E(Ecur, idx[0], idx[1], x, y)
if (E2 < Ebest):
Ecur, x[idx] = E2, new
Ebest = E2  # update Ebest
else:
Ecur, x[idx] = E1, old

# record time and Ebest of this iteration
end_time = time.time()
energy_record[0].append(end_time - initial_time)
energy_record[1].append(Ebest)

return x, energy_record

def sign(data, translate):
"""Map a dictionary for the element of data.

Example:
To convert every element in data with value 0 to -1, 255 to 1,
use signed = sign(data, {0: -1, 255: 1})
"""
temp = np.array(data)
return np.vectorize(lambda x: translate[x])(temp)

In [21]: beta, eta, h = 1e-3, 2.1e-3, 0.0
In [22]: im = Image.open('flipped.png')
In [23]: data = sign(im.getdata(), {0: -1, 255: 1})  # convert to {-1, 1}
In [24]: y = data.reshape(im.size[::-1]) # to matrix
In [25]: E, localized_E = E_generator(beta, eta, h)
In [26]: result, energy_record = ICM(y, E, localized_E)
In [27]: result = sign(result, {-1: 0, 1: 255})
In [28]: output_image = Image.fromarray(result).convert('1', dither=Image.NONE)

output_image.show()

In [30]: plt.plot(*energy_record)
Out[30]: [<matplotlib.lines.Line2D at 0xa63ecd0>]

In [31]: plt.xlabel('Time(s)')
Out[31]: <matplotlib.text.Text at 0xa5efb70>

In [32]: plt.ylabel('Energy')
Out[32]: <matplotlib.text.Text at 0xa5f3b70>

In [33]: plt.show()

In [60]: original = Image.open('in.png')
In [61]: x = np.reshape(original.getdata(), original.size[::-1])
In [62]: 1 - np.count_nonzero(x - result) / float(x.size)
Out[62]: 0.9621064814814815

### 模拟退火

def temperature(k, kmax):
"""Schedule the temperature for simulated annealing."""
return 1.0 / 500 * (1.0 / k - 1.0 / kmax)

def prob(E1, E2, t):
"""Probability transition function for simulated annealing."""
return 1 if E1 > E2 else np.exp((E1 - E2) / t)

def simulated_annealing(y, kmax, E, localized_E, temp_dir):
"""Simulated annealing process for image denoising.

Parameters
----------
y: array_like
The noisy binary image matrix ranging in {-1, 1}.
kmax: int
The maximun number of iterations.
E: function
Energy function.
localized_E: function
Localized version of E.
temp_dir: path
Directory to save temporary results.

Returns
----------
x: array_like
The denoised binary image matrix ranging in {-1, 1}.
energy_record:
[time, Ebest] records for plotting.
"""
x = np.array(y)
Ebest = Ecur = E(x, y)  # initial energy
initial_time = time.time()
energy_record = [[0.0, ], [Ebest, ]]

for k in range(1, kmax + 1):  # iterate kmax times
start_time = time.time()
t = temperature(k, kmax + 1)
print "k = %d, Temperature = %.4e" % (k, t)
accept, reject = 0, 0  # record the acceptance of alteration
for idx in np.ndindex(y.shape):  # for each pixel in the matrix
old, new, E1, E2 = localized_E(Ecur, idx[0], idx[1], x, y)
p, q = prob(E1, E2, t), random()
if p > q:
accept += 1
Ecur, x[idx] = E2, new
if (E2 < Ebest):
Ebest = E2  # update Ebest
else:
reject += 1
Ecur, x[idx] = E1, old

# record time and Ebest of this iteration
end_time = time.time()
energy_record[0].append(end_time - initial_time)
energy_record[1].append(Ebest)

print "k = %d, accept = %d, reject = %d" % (k, accept, reject)
print "k = %d, %.1f seconds" % (k, end_time - start_time)

# save temporary results
temp = sign(x, {-1: 0, 1: 255})
temp_path = os.path.join(temp_dir, 'temp-%d.png' % (k))
Image.fromarray(temp).convert('1', dither=Image.NONE).save(temp_path)
print "[Saved]", temp_path

return x, energy_record

time-energy 关系图

### 后话

posted @ 2015-02-01 00:42  Joyee  阅读(...)  评论(...编辑  收藏