画像のノイズリダクション¶

このサンプルコードでは、Fixstars Amplify を用いて、ノイズが加わった以下の画像から、ノイズが加わる前の元の画像を推定します。

定式化の概要¶

今回は簡単のため、白黒画像を取り扱います。ノイズ除去後の画像とノイズの加わった画像の関係として、以下を仮定します。

ノイズ除去後の画像とノイズの入った画像は一致することが多い
ノイズ除去後の画像では隣り合う画素は同じ色であることが多い

これらの条件をできるだけ満たすように最適化を行うことで、ノイズの入った画像からノイズを除去します。

定数と変数¶

白黒画像の各画素は黒または白のいずれかなので、画像の各画素を $-1$ または $1$ の値で表現することができます。以下では、黒は $-1$ に対応し、白は $1$ に対応するということにします。

ノイズが加わった画像について、画素の色を表す -1 または 1 の数値を $y_{ij}$ とおきます。

$$ y_{ij} = \begin{cases} -1 \quad \text{(対応する画素が黒色)} \\ +1 \quad \text{(対応する画素が白色)} \end{cases} $$

ノイズ除去後の画像については、各画素の値があらかじめ分かっているわけではないですが、黒または白の値を取ることは確実なので、イジング変数 ($\pm1$ の値を取る変数) $s_{ij}$ で表します。このイジング変数について最適化を行うことで、ノイズの除去を行うことができます。

$$ s_{ij} \in \{-1, 1\} $$

目的関数¶

ノイズ除去後の画像に期待される条件は、以下の 2 種類です。

ノイズ除去後の画像とノイズの入った画像をできるだけ一致させたい
ノイズ除去後の画像において、隣り合う画素はできるだけ同じ色にしたい

これらの条件を多項式の最小化問題として書き下します。

まず、「ノイズ除去後の画像とノイズの入った画像をできるだけ一致させたい」という条件を考えます。ノイズ除去後の画像とノイズの入った画像の 2 つの画像において、同じ位置 (i, j) にある画素が一致する場合 $-y_{ij}s_{ij} = -1$、一致しない場合 $-y_{ij}s_{ij} = 1$ となるので、この式を画素ごとに足し合わせた多項式

$$ f_1 = \sum_{i, j} -y_{ij} s_{ij} $$

の値が小さければ小さいほど、ノイズ除去後の画像とノイズの入った画像の一致度が高くなるといえます。

次に、「ノイズ除去後の画像において、隣り合う画素はできるだけ同じ色にしたい」という条件を考えます。隣り合う位置にある 2 つの画素 $s_{ij}$, $s_{i'j'}$ について、これらの色が一致する場合 $-s_{ij}s_{i'j'} = -1$、一致しない場合 $-s_{ij}s_{i'j'} = 1$ となるので、この式を隣り合う画素のペアごとに足し合わせた多項式

$$ f_2 = \sum_{s_{i, j} \text{と} s_{i', j'} \text{が隣り合う}} -s_{i, j} s_{i', j'} $$

の値が小さければ小さいほど、隣り合う画素が同じ色になっている割合が高くなるといえます。

これらの 2 式 $f_1$ と $f_2$ を重みをつけて足し合わせることで、2 つの条件を両方ともできるだけ満たす場合に小さな値をとるような目的関数 $f$ を構築できます。

$$ \begin{align} f & = f_1 + \eta f_2\\ f_1 &= \sum_{i, j} -y_{ij} s_{ij}, \\ f_2 &= \sum_{s_{i, j} \text{と} s_{i', j'} \text{が隣り合う}} -s_{i, j} s_{i', j'} \end{align} $$

ここで、$\eta>0$ というパラメータを導入しました。これにより $f_1$ と $f_2$ の強さの調整が出来ます。$\eta$ を小さくすると 2 つの条件のうち「ノイズ除去後の画像とノイズの入った画像をできるだけ一致させたい」割合が強くなり、大きくすると「隣り合う画素をできるだけ同じ色にしたい」割合が強くなります。

参考¶

ノイズが乗った画像の読み込み¶

以下のようなノイズの乗った白黒画像 (noisy.png) を用意しました。

ノイズの乗った画像

この画像に対してノイズ除去を行い、元の画像を推定することが今回の目的となっています。

まず、ノイズが乗った画像を読み込み、NumPy 配列に変換します。

In [ ]:

from PIL import Image
import numpy as np
    
# 元画像のイジング配列を作成
img = Image.open("noisy.png")
noisy_img = np.where(np.array(img) >= 128, 1, -1)

NumPy 配列 noisy_img は画像サイズと同じく 81 x 196 の二次元配列であり、各要素は -1 (黒) または 1 (白) となっています。この NumPy 配列を Matplotlib を用いて可視化すると以下のようになります。

In [ ]:

%matplotlib inline

In [ ]:

import matplotlib.pyplot as plt
    
print(noisy_img.shape)
plt.imshow(noisy_img, cmap="gray")
plt.show()

イジング変数配列の作成¶

次に、ノイズ除去後の画像の各画素の値を表すイジング変数の配列 s を生成します。s はノイズの乗った画像と同じサイズ、つまり NumPy 配列 noisy_img と同じサイズの二次元配列です。

変数の生成には VariableGenerator を使います。VariableGenerator の array メソッドを用いると、多次元配列の形で変数配列を作成できるので、目的関数の計算に便利です。

In [ ]:

from amplify import VariableGenerator
    
gen = VariableGenerator()
# 元の画像と同じ 2 次元配列の形にイジング変数配列を生成
s = gen.array("Ising", shape=noisy_img.shape)

目的関数¶

ノイズの乗った画像に対応する NumPy 配列 noisy_img と元の画像を表すイジング変数配列 $s$ を用いて、目的関数を構築します。

まず、「ノイズ除去後の画像とノイズの入った画像は一致することが多い」という条件を表す関数を作成します。この関数は $ f_1 = \sum_{i, j} -y_{ij} s_{ij} $ で表され、一致する画素の数が多ければ大きいほど小さい値を取ります。

In [ ]:

f1 = -(noisy_img * s).sum()

次に、「ノイズ除去後の画像において、隣り合う画素はできるだけ同じ色にしたい」という条件を表す関数を作成します。この関数は以下の式で表され、同じ色の隣り合う画素が多ければ多いほど小さな値を取ります。

$$ f_2 = \sum_{s_{i, j} \text{と} s_{i', j'} \text{が隣り合う}} -s_{i, j} s_{i', j'} $$

In [ ]:

f2 = -((s[:, :-1] * s[:, 1:]).sum() + (s[:-1, :] * s[1:, :]).sum())

ここで、上の f2 の式の第一項は左右に隣り合う画素のペアすべてについて、対応するイジング変数の積 $s_{i, j}s_{i', j'}$ の総和を取ったものとなっています。s[:, :-1] と s[:, 1:] は同じ形を持つ二次元配列であり、これら 2 つの配列において同じ位置にある 2 つのイジング変数は s において左右に隣り合っていることに注意してください。同様に、第二項は上下に隣り合う画素のペアすべてについて対応するイジング変数の積の総和を取ったものです。

作成した 2 つの関数を適切な重みを付けて足し合わせ、目的関数を作成します。今回は第一項 f1 に対する第二項 f2 の重みを 0.4 とします。この目的関数ができるだけ小さくなるように変数配列 s の値を決定したとき、s が表す画像はノイズが加わる前の画像に近くなることが期待されます。

In [ ]:

# 強度パラメータ
eta = 0.4
    
# 上記 2 つの関数の値が同時に小さくなるように目的関数を設定する
objective = f1 + eta * f2

ソルバークライアントの設定とマシンの実行¶

組合せ最適化ソルバーとして、Amplify AE を使用します。Amplify AE に対応するソルバークライアント (FixstarsClient) を作成し、パラメータを設定します。

In [ ]:

from amplify import FixstarsClient
from datetime import timedelta
    
# クライアントの設定
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=2000)
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境で使用する場合は、Amplify AEのアクセストークンを入力してください

作成した目的関数とソルバークライアントを用いて、組合せ最適化を実行します。以下のセルを実行すると、Amplify AE による目的関数の最小化が行われます。

In [ ]:

from amplify import solve
    
result = solve(objective, client)

解の取得と結果の表示¶

得られたうち最良の解における変数配列 s の値を取得します。

In [ ]:

# イジング変数に解を代入
s_values = s.evaluate(result.best.values)

s_values は s と同じサイズの NumPy 2 次元配列であり、ノイズ除去後の画像の各画素の値を表しています。

最後に、s_values が表す画像を表示してみましょう。最初の画像 (noisy.png) と比較するとノイズが減少したことが確認できます。

In [ ]:

plt.imshow(s_values, cmap="gray")  # 復元画像
plt.show()