目次
画像処理におけるスパースモデリング
画像のスパースモデリングとは スパース：「疎ら」 なモデリングを行う，ということです．このようなスパースモデリングは，画像処理にとどまらずノイズ除去や，領域のセグメンテーションなどに用いられています．
ここでは，画像処理において，スパースモデリングとはどういうことを意味するのか，pythonを用いてコーディングを行います．
1. L1 正則化に基づく画像の再構成
L1画像の再構成では，具体的に以下のような評価関数の最小化を考えます． Input の画像から Outputの画像に変換します．
\begin{align} F = \frac{1}{2} |I_{input} - I_{output}|^2 + \lambda |I_{output}|_{1} \end{align}
この評価関数の第一項目は入力画像と出力画像の2乗誤差です．これは，Input画像とOutput画像が明らかに違うようなものは求めてない，ということです．
大事なのは第二項目です．この項は$I_{input}$のL1ノルムに対応します．L1ノルムの最小化はスパースな解，つまり要素として0を多く含む解が得られるようになります．今回の場合，$I_{output}$のL1ノルムの最小化を行っているので，$I_{output}$の画素の多くが0になるような Output を得ることができます．
なお，0ではなく別の値$A$の値をスパースにしたいときなどは，上の評価関数を
\begin{align} F = \frac{1}{2} |I_{input} - I_{output}|^2 + \lambda |I_{output} - A|_{1} \end{align} とすることで解決できます．上のように，バックグラウンドの値は白なので，GrayScale の画像では$A=255$とする必要がありますね．
2. Total Variationに基づく画像の再構成
上のL1ノルムの最小化は，バックグラウンドのノイズを除去するのに使われます．これは，背景部分がほとんどある値$A$になっているだろう，という仮定によるものでした．
それでは次に，以下のような評価関数を考えてみましょう．
\begin{align} F’ = \frac{1}{2} |I_{input} - I_{output}|^2 + \lambda |\nabla I_{output}| \end{align}
ただし， $ |\nabla I_{output}| = | \nabla_x I_{output}| + | \nabla_y I_{output}| $ の和であり，$ \nabla_x I_{output} $は$ I_{output} $の隣接するピクセルの$x$方向の差分です．この正則化をTotal Variation (以下TV)正則化と呼びます．
今度の正則化はOutput 画像の画素の隣接ピクセルの差分に関する正則化です．差分のL1をとっているため，差分がスパースつまり画素値の変化が画像内で少なくなるような効果を生み出します．
具体的には以下のような画像のノイズを削減しようとするときに，この正則化を用いることで解決できます．
Python による実装
それでは，この最適化をプログラムしてみましょう． ただし，L1ノルム正則化の最適化問題は，簡単には実装できません． 最適化問題を解く方法として，様々な方法が提案されていますが，この記事では Split Bregman 法 により最適化を行います．
1. L1正則化の場合
評価関数を再度記述すると以下のようになります．
\begin{align} F = \frac{1}{2} |I_{input} - I_{output}|^2 + \lambda |I_{output}|_{1} \end{align}
これを最小化する解は，実は陽に書くことができます．それを具体的に示していきましょう．
評価関数の変数を単純化します． \begin{align} F = \frac{1}{2} |x - b|^2 + \lambda |x|_{1} \end{align}
上の評価関数を要素ごとの和で書き表すと
\begin{align} F = \frac{1}{2}\sum_{i} (x_i - b_i)^2 + \lambda \sum_{i} |x_{i}| \end{align}
が得られます．よって，各変数$x_i$に関して最小化を行えばよいことが分かります．
次に，$x_i=a_i e_i$，$a_i>0,\ e_i=1,\ or, \ -1$とします．つまり，$x_i$を符号のIndicator $e_i$と大きさを表す$a_i$に分割します．これを用いると，上の式は
\begin{align} F = \frac{1}{2}\sum_{i} (a_i e_i - b_i)^2 + \lambda \sum_{i} a_i \end{align} となります．
ここで，$a_i>0$に注意します．$(a_i e_i - b_i)^2$をできるだけ小さくしたいのですが，$b_i$と$e_i$とがもし異符号であれば，$(a_ie_i - b_i)^2 > b_i^2$となります．よって，最小化したいのであれば，$e_i = {\rm sgn}(b_i)$となることが分かります．
このように$e_i$が求まれば，あとは$a_i$に関する最小化を行うことで解決します．具体的には
\begin{align} F = \frac{1}{2}\sum_{i} (a_i {\rm sgn}(b_i) - b_i)^2 + \lambda \sum_{i} a_i \end{align} \begin{align} = \frac{1}{2}\sum_{i} (a_i - |b_i|)^2 + \lambda \sum_{i} a_i \end{align} \begin{align} = \frac{1}{2}\sum_{i} (a_i - |b_i| + \lambda )^2 + 2\lambda|b_i| - \lambda^2 \end{align}
これから，$a_i>0$に注意して，
\begin{align} a_i = {\rm max} (|b_i|-\lambda,0) \end{align} という解が得られます．
つまり．L1ノルムの制約付き最小化問題では，このように設定することで，解が得られることが分かります．正則化パラメータを調整することで，L1ノルムの最適化問題を解くことができます．
実装
ImageJのサンプル画像から拾ってきました．
この画像にノイズを加えたデータを用いて，ノイズ除去を行ってみます．
import numpy as np import cv2 import matplotlib.pyplot as plt import random
img_path = "Blood2.jpg"
img_load = cv2.imread(img_path)
I_t = cv2.cvtColor(img_load, cv2.COLOR_RGB2GRAY)## add noise mu = np.mean(I_t) sigma = np.std(I_t) dB = 10 I_noise = 10**(-dB/20)*np.reshape([random.gauss(mu, sigma) for i in range(np.size(I_t))], np.shape(I_t)) I = I_t + I_noise max_I = np.max(I) min_I = np.min(I) I = np.round((I - min_I)*255/(max_I - min_I))
plt.subplot(1, 2, 1) plt.gray() plt.imshow(I_t) plt.subplot(1,2,2) plt.gray() plt.imshow(I)
Noisy データの完成
LAMBDA = 30 ## 正則化パラメータ THRE = 200 ## だいたい画素値200 くらいが background B = I - THRE t = np.sign(B)*B - LAMBDA I_reconst = t*(t>0)*np.sign(B) + THRE
plt.subplot(1, 3, 1) plt.gray() plt.imshow(I_t) plt.axis("off") plt.title("Original") plt.subplot(1,3,2) plt.gray() plt.imshow(I) plt.axis("off") plt.title("Noisy") plt.subplot(1,3,3) plt.gray() plt.imshow(I_reconst) plt.axis("off") plt.title("Reconstructed")
それでは，正則化パラメータを変えてみましょう．
LAMBDA_LIST = [5, 30,50,100] THRE = 200 I_LIST = [] for Lambda in LAMBDA_LIST: B = I - THRE t = np.sign(B)*B - Lambda I_reconst = t*(t>0)*np.sign(B) + LAMBDA I_LIST.append(I_reconst)
for i in range(len(LAMBDA_LIST)): plt.subplot(2, len(LAMBDA_LIST)/2, i+1) plt.gray() plt.imshow(I_LIST[i]) plt.axis("off") plt.title("$\lambda = $" + str(LAMBDA_LIST[i]))
正則化パラメータの決定方法はT-curve 法やL-curve 法，Cross Varidation に基づく方法などがあります．またこれに関しても紹介していけたらと思っています．
長くなったので，今回の記事はここで一旦終了します． 次の記事からは，微分のL1ノルム最小化，つまりTotal Variation の正則化項を加えた場合のイメージングアプローチの最小化を行っていきます．次の記事は下です．
スパースモデリングに基づく画像の再構成 Part2. Total Variation最小化(Split Bregman)に基づく画像再構成 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-|LPixel(エルピクセル)
参考文献
https://imagej.nih.gov/ij/
[2] THE SPLIT BREGMAN METHOD FOR L1 REGULARIZED PROBLEMS
(TOM GOLDSTEIN, STANLEY OSHER)
ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf