当社はCookieを使用して、お客様が当社のWebサイトでより良い体験を得られるようにしています。引き続き閲覧する場合は、当社グループのプライバシーポリシー 
に同意したことになります。
Jupyter Notebookで実行する
サンプルコードアニーリングマシンを用いてノイズが加わった画像から元の画像の推定を行うことを考えます。
下記の仮定に基づいてノイズ除去を試みます。
ここでは簡単のため白黒画像を取り扱います。
画素のデータは黒と白の二値での表現が可能なので、各画素の値を二値変数を用いて表すことができます。
上記の仮定を表現する画素同士の相互作用を目的関数を用いて表現し、これを最適化することによって元の画像を推定することができます。
画素の集合を $V$ とし、各画素を表すインデックスを $i\in V$ とします。
まず、入力画素を表すイジング変数を $y$ としそれぞれの色に対応する値を以下のように表します。
また、出力画素に対応した二値のイジング変数を以下のように表します。
$$ s_{i} = \left\{ \begin{align} &+1 \quad\text{(白)}\\ &-1 \quad \text{(黒)} \end{align} \right. \quad i\in V\\ $$入力画像と出力画像は概ね一致するという仮定 (ノイズがそれほど多くないという仮定) により、入力画素と出力画素は同じ値になるようにします。つまり、$s_i$ と $y_i$ は同じ値を持つときに値が小さくなるような目的関数を導入します。例えば以下のように与えられます。
$$ f_1 = - \sum_{i\in V} y_{i} s_{i} $$$y_{i}$ と $s_{i}$ が同じ値を持つと 上記の目的関数の値は減少し、異なった値を持つと増加するので、全ての $i\in V$ において $y_{i} = s_{i}$ である場合に $f_1$ は最小値をとります。しかしながら、入力画像にはノイズがのっているので、出力画像が入力画像と同じになってしまうとノイズを減らすことができません。
そこで、隣り合う画素は同じ色になりやすいという仮定を考慮します。
つまり、隣り合った出力画素が同じ値を持つ場合に値が小さくなるような目的関数を導入します。例えば以下のように与えられます。
ここで隣接する画素のペアの集合を $E$ としました。全ての出力画素が同じ値を持つと $f_2$ は最小の値をとります。しかし、全ての画素が同じ値になってしまうと全てが白または黒の画像になってしまうので、元の画像の情報が失われてしまいます。
そこで、$f_1$ と $f_2$ を適切に足し合わせることで、出力画像が入力画像と近い値をとりつつノイズと思われる画素のみ除去することを試みます。
$$ \begin{align} f & = f_1 + \eta f_2\\ &=- \sum_{i\in V}y_is_i - \eta \sum_{(i,j)\in E}s_i s_j \end{align} $$ここで、$\eta>0$ というパラメータを導入しました。これにより $f_1$ と $f_2$ の強さの調整が出来ます。$\eta$ が大きいほどノイズ除去を行う項が強いことを意味しています。
この目的関数を最小化しイジング変数 $s$ の値を画素の値と解釈することで、ノイズを除去した画像が得られます。
まずは、画像データをダウンロードする関数と、ダウンロードした画像をイジング変数配列に変換する関数を定義します。
import os
from PIL import Image
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
# 元画像のイジング配列を作成
img = Image.open("sample.png")
x = np.where(np.array(img) >= 128, 1, -1)
plt.imshow(x, cmap="gray")
plt.show()
print(x.shape)
print(x)
次に、画素を無作為に選びその値を反転することでノイズを表現する関数を定義します。
# 画像にノイズを追加する関数
def get_noisy_img_array(img_array):
# 2次元配列を1次元に変換して扱いやすくする
img_shape = img_array.shape
flattened_img = img_array.flatten()
# 最大値と最小値を入れ替える関数を定義
min_v = min(flattened_img)
max_v = max(flattened_img)
def invert_value(v):
return min_v + max_v - v
# ノイズの割合
ratio = 0.02
# ノイズをのせる画素をランダムに選択して反転
for idx in np.random.choice(len(flattened_img), int(ratio * len(flattened_img))):
flattened_img[idx] = invert_value(flattened_img[idx])
# 元の配列の形に戻す
return flattened_img.reshape(*img_shape)
# ノイズ画像のイジング配列を作成
y = get_noisy_img_array(x)
plt.imshow(y, cmap="gray")
plt.show()
print(y.shape)
print(y)
次に、イジング変数の配列 s を生成します。入力画像のデータ y を $h \times w$ の2次元配列とすると、出力画像に対応するイジング変数
s も 同じく $h \times w$ の2次元配列となります。
変数の生成には IsingSymbolGenerator を使います。IsingSymbolGenerator の array
メソッドでは、入力画像のデータ y と同じ配列の形で変数配列を作成できるので、目的関数の計算に便利です。
from amplify import IsingSymbolGenerator, IsingPoly, sum_poly
# 画像の高さ(h),幅(w)を取得
h, w = y.shape
gen = IsingSymbolGenerator()
s = gen.array(h, w) # h x w の配列の形にイジング変数を生成
入力画像データの配列 $y$ と出力画像に対応したイジング変数配列 $s$ を用いて、目的関数を構築します。
# 強度パラメータ
eta = 0.333
# 目的関数 f を計算
# - \sum_{i\in V} y_{i} s_{i}
f1 = sum_poly(-s * y)
# -\sum_{(i,j)\in E} s_i s_j
f2 = sum_poly(
h - 1, lambda i: sum_poly(w, lambda j: -s[i, j] * s[i + 1, j])
) + sum_poly(h, lambda i: sum_poly(w - 1, lambda j: -s[i, j] * s[i, j + 1]))
f = f1 + eta * f2
次にクライアントを設定し、先ほど与えた目的関数の最小値に対応する解をイジングマシンで探索します。
from amplify.client import FixstarsClient
from amplify import Solver
# クライアントの設定
client = FixstarsClient()
client.parameters.timeout = 1000 # タイムアウト1秒
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # ローカル環境で使用する場合は、Amplify AEのアクセストークンを入力してください
# ソルバの設定と結果の取得
solver = Solver(client)
result = solver.solve(f)
# 解の取得
values = result[0].values
# イジング変数に解を代入
output = s.decode(values, 1)
plt.imshow(output, cmap="gray") # 復元画像
plt.show()
plt.imshow(x, cmap="gray") # 元画像
plt.show()
plt.imshow(y, cmap="gray") # ノイズ画像
plt.show()
